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 tesselating elements in such a way that each integration point is represented by one subsection of a cell.

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 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 ip_tesselated()

from pathlib import Path
from tempfile import mkdtemp

import ogstools as ogs
from ogstools import examples
from ogstools.meshlib.gmsh_meshing import rect
from ogstools.msh2vtu import msh2vtu

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

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

tmp_dir = Path(mkdtemp())
mesh_path = tmp_dir / "mesh.msh"
vtu_path = tmp_dir / "mesh_domain.vtu"


def simulate_and_plot(elem_order: int, quads: bool, intpt_order: int):
    rect(
        lengths=1,
        n_edge_cells=6,
        structured_grid=quads,
        order=elem_order,
        out_name=mesh_path,
    )
    msh2vtu(mesh_path, tmp_dir, log_level="ERROR")

    model = ogs.Project(
        output_file=tmp_dir / "default.prj",
        input_file=examples.prj_mechanics,
    )
    model.replace_text(intpt_order, xpath=".//integration_order")
    model.write_input()
    model.run_model(write_logs=True, args=f"-m {tmp_dir} -o {tmp_dir}")
    mesh = ogs.MeshSeries(tmp_dir / "mesh.pvd").mesh(-1)
    int_pts = mesh.to_ip_point_cloud()
    ip_mesh = mesh.to_ip_mesh()

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

Triangles with increasing integration point order

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)

• 
• 

OGS finished with project file /tmp/tmpd2zf456y/default.prj.
Execution took 0.1532599925994873 s
Project file written to output.
/builds/ogs/tools/ogstools/.venv-devcontainer/lib/python3.10/site-packages/pyvista/core/pointset.py:843: PyVistaDeprecationWarning: `PolyData` constructor parameter `n_faces` is deprecated and no longer used.
  warnings.warn(

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

• 
• 

OGS finished with project file /tmp/tmpd2zf456y/default.prj.
Execution took 0.15220999717712402 s
Project file written to output.
/builds/ogs/tools/ogstools/.venv-devcontainer/lib/python3.10/site-packages/pyvista/core/pointset.py:843: PyVistaDeprecationWarning: `PolyData` constructor parameter `n_faces` is deprecated and no longer used.
  warnings.warn(

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

• 
• 

OGS finished with project file /tmp/tmpd2zf456y/default.prj.
Execution took 0.1540071964263916 s
Project file written to output.
/builds/ogs/tools/ogstools/.venv-devcontainer/lib/python3.10/site-packages/pyvista/core/pointset.py:843: PyVistaDeprecationWarning: `PolyData` constructor parameter `n_faces` is deprecated and no longer used.
  warnings.warn(

Quadratic triangles

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

• 
• 

OGS finished with project file /tmp/tmpd2zf456y/default.prj.
Execution took 0.16376781463623047 s
Project file written to output.
/builds/ogs/tools/ogstools/.venv-devcontainer/lib/python3.10/site-packages/pyvista/core/pointset.py:843: PyVistaDeprecationWarning: `PolyData` constructor parameter `n_faces` is deprecated and no longer used.
  warnings.warn(

Quadrilaterals with increasing integration point order

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)

• 
• 

OGS finished with project file /tmp/tmpd2zf456y/default.prj.
Execution took 0.1508471965789795 s
Project file written to output.
/builds/ogs/tools/ogstools/.venv-devcontainer/lib/python3.10/site-packages/pyvista/core/pointset.py:843: PyVistaDeprecationWarning: `PolyData` constructor parameter `n_faces` is deprecated and no longer used.
  warnings.warn(

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

• 
• 

OGS finished with project file /tmp/tmpd2zf456y/default.prj.
Execution took 0.15270280838012695 s
Project file written to output.
/builds/ogs/tools/ogstools/.venv-devcontainer/lib/python3.10/site-packages/pyvista/core/pointset.py:843: PyVistaDeprecationWarning: `PolyData` constructor parameter `n_faces` is deprecated and no longer used.
  warnings.warn(

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

• 
• 

OGS finished with project file /tmp/tmpd2zf456y/default.prj.
Execution took 0.15684008598327637 s
Project file written to output.
/builds/ogs/tools/ogstools/.venv-devcontainer/lib/python3.10/site-packages/pyvista/core/pointset.py:843: PyVistaDeprecationWarning: `PolyData` constructor parameter `n_faces` is deprecated and no longer used.
  warnings.warn(

Quadratic quadrilateral

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

• 
• 

OGS finished with project file /tmp/tmpd2zf456y/default.prj.
Execution took 0.1574690341949463 s
Project file written to output.
/builds/ogs/tools/ogstools/.venv-devcontainer/lib/python3.10/site-packages/pyvista/core/pointset.py:843: PyVistaDeprecationWarning: `PolyData` constructor parameter `n_faces` is deprecated and no longer used.
  warnings.warn(