"""
Analyzing integration point data
================================

This examples shall demonstrate how we can better visualize integration point
data (raw data used in OGS's equation system assembly without output related
artefacts), by tessellating elements in such a way that each integration point is
represented by one subsection of a cell.
"""

# %% [markdown]
# For brevity of this example we wrap the entire workflow from meshing and
# simulation to plotting in a parameterized function. The main function of
# importance is :func:`~ogstools.mesh.to_ip_mesh`. We will use this
# function to show the tessellated visualization of the integration point data
# for different element and integration point orders and element types. Note,
# you can also tessellate an entire MeshSeries via
# :meth:`~ogstools.MeshSeries.ip_tesselated`


import ogstools as ot
from ogstools import examples

ot.plot.setup.dpi = 75
ot.plot.setup.show_element_edges = True

sigma_ip = ot.variables.stress.replace(
    data_name="sigma_ip", output_name="IP_stress"
)


def simulate_and_plot(elem_order: int, quads: bool, intpt_order: int):

    meshes = ot.Meshes.from_gmsh(
        ot.gmsh_tools.rect(
            lengths=1,
            n_edge_cells=6,
            structured_grid=quads,
            order=elem_order,
        ),
        log=False,
    )

    prj = ot.Project(
        input_file=examples.prj_mechanics,
    ).copy()
    prj.replace_text(intpt_order, xpath=".//integration_order")
    model = ot.Model(prj, meshes)
    sim = model.run()
    mesh = sim.meshseries.mesh(-1)
    int_pts = ot.mesh.to_ip_point_cloud(mesh)
    ip_mesh = ot.mesh.to_ip_mesh(mesh)

    fig = ot.plot.contourf(mesh, ot.variables.stress)
    fig.axes[0].scatter(
        int_pts.points[:, 0], int_pts.points[:, 1], color="k", s=10
    )
    fig = ot.plot.contourf(ip_mesh, sigma_ip)
    fig.axes[0].scatter(
        int_pts.points[:, 0], int_pts.points[:, 1], color="k", s=10
    )


# %% [markdown]
# Triangles with increasing integration point order
# -------------------------------------------------
# .. dropdown:: Why does the stress not change with the integration point order?
#
#     In linear triangular finite elements, the shape functions used
#     to interpolate displacements are linear functions of the coordinates.
#     As this is a linear elastic example, the displacements are linear.
#     The strain, which is obtained by differentiating the displacement, will
#     thus be constant throughout the element. The stress, which is related to
#     the strain through a constitutive relationship will also be constant
#     throughout the element. Thus, the stress is not affected by the
#     integration point order.

simulate_and_plot(elem_order=1, quads=False, intpt_order=2)

# %%
simulate_and_plot(elem_order=1, quads=False, intpt_order=3)

# %%
simulate_and_plot(elem_order=1, quads=False, intpt_order=4)

# %% [markdown]
# Quadratic triangles
# -------------------

simulate_and_plot(elem_order=2, quads=False, intpt_order=4)

# %% [markdown]
# Quadrilaterals with increasing integration point order
# ------------------------------------------------------
# .. dropdown:: Why does the stress change here?
#
#     In contrast to triangular elements, quadrilateral elements use bilinear
#     shape functions. Thus, the differentiation of the displacement leads to
#     bilinear strain. The stress in turn is bilinear as well and can change
#     within the elements. The number of integration points consequently
#     affects the resulting stress field.

simulate_and_plot(elem_order=1, quads=True, intpt_order=2)

# %%
simulate_and_plot(elem_order=1, quads=True, intpt_order=3)

# %%
simulate_and_plot(elem_order=1, quads=True, intpt_order=4)

# %% [markdown]
# Quadratic quadrilateral
# -----------------------

simulate_and_plot(elem_order=2, quads=True, intpt_order=4)