Mesh Tutorial

This page is based on a Jupyter notebook.

Goals

In this tutorial the user will learn to use some of the tools provided with OpenGeoSys to create simple mesh for a simulation.

Requirements

A running version of OpenGeoSys (it can be installed via pip, as a container or as ‘CLI with utilities’ binary download)
It is assumed, that regardless of how OGS is installed, it is added to the environment variables and is accessible from the shell.
In case of a problem, please consult the installation and troubleshooting pages.
For viewing a mesh with PyVista, Python is required.
To run the Jupyter notebook that accompanies this tutorial, a Python environment with Jupyter, PyVista, NumPy and Matplotlib packages is needed.

Please note that this tutorial doesn’t cover Python. The user is expected to have basic understanding of it. Alternatively, other tools capable of reading mesh files can be used (for example: ParaView).

A basic understanding of meshes would make following this tutorial easier but is not required.

Conventions

If parameters can take different values (for example to apply them in different directions) variable part of command will be indicated by $ sign followed by name.

For example directional parameters set -dx -dy -dz will be indicated by -d$direction $value.
This doesn’t apply to code snippets, there the parameters will be written explicitly.

! is used to run shell commands from within the Jupyter notebook.

Useful links

It is not mandatory to read the following articles before starting with the tutorial, but they can provide some expansion and additional explanation to the content presented here.

Overview over the tutorial

In this tutorial, some tools that can be used to create a mesh will be presented with a short description and an example. These tools are shipped as executables with OGS:

Afterwards, you can test your understanding with a couple of exercises.

Preparations

We first set an output directory and import some required libraries. This tutorial uses PyVista for visualization of the meshes.

import os… (click to toggle)

import os
from pathlib import Path

# All generated files will be put there
out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
    out_dir.mkdir(parents=True)

# Use pyvista for visualization of the generated meshes… (click to toggle)

# Use pyvista for visualization of the generated meshes
import pyvista as pv
pv.set_plot_theme("document")
pv.set_jupyter_backend("static")

import matplotlib.pyplot as plt
import numpy as np

# Define some plotting functions… (click to toggle)

# Define some plotting functions

def show_mesh(mesh):
    plotter = pv.Plotter()
    plotter.add_mesh(mesh, show_edges=True)
    plotter.view_xy()
    plotter.add_axes()
    plotter.show_bounds(mesh, xtitle="x", ytitle="y")
    plotter.window_size = [800, 400]
    plotter.show()

def show_mesh_3D(mesh):
    plotter = pv.Plotter()
    plotter.add_mesh(mesh, show_edges=True)
    plotter.add_axes()
    plotter.window_size = [800, 400]
    plotter.show()

def show_cell_sizes_x(mesh):
    line = mesh.slice_along_line(pv.Line((0,1e-8,0),(7,1e-8,0), resolution=2))

    fig, ax = plt.subplots()

    cell_lengths_x = np.diff(line.points[:,0])
    num_cells_x = len(cell_lengths_x)
    cell_idcs_x = np.arange(1, num_cells_x+1)

    ax.plot(cell_idcs_x, cell_lengths_x, marker=".", label="cell sizes from mesh")
    ax.plot(cell_idcs_x, 1.2**cell_idcs_x * 0.5, label="cell sizes expected, up to a factor")

    ax.legend()
    ax.set_xtitle("cell id in $x$ direction")
    ax.set_ytitle("cell size in $x$ direction")

    ax.set_yscale("log")

Structured mesh generation

You can start by checking the version of the generateStructuredMesh tool, it will test whether your OpenGeoSys is correctly installed and available:

! generateStructuredMesh --version
assert _exit_code == 0

generateStructuredMesh  version: 6.5.3-407-g9a1da6a18

This command should deliver the following output:

generateStructuredMesh  version: - some version -

You can view basic documentation using this command:

! generateStructuredMesh --help
assert _exit_code == 0

USAGE: 

   generateStructuredMesh  [-h] [--version] [--] -e <line|tri|quad|hex
                           |prism|tet|pyramid> -o <file name of output
                           mesh> [--lx <real>] [--ly <real>] [--lz <real>]
                           [--nx <integer>] [--ny <integer>] [--nz
                           <integer>] [--dx0 <real>] [--dy0 <real>] [--dz0
                           <real>] [--dx-max <real>] [--dy-max <real>]
                           [--dz-max <real>] [--mx <real>] [--my <real>]
                           [--mz <real>] [--ox <real>] [--oy <real>] [--oz
                           <real>]


Where: 

   -h,  --help
     Displays usage information and exits.

   --version
     Displays version information and exits.

   --,  --ignore_rest
     Ignores the rest of the labeled arguments following this flag.

   -e <line|tri|quad|hex|prism|tet|pyramid>,  --element-type <line|tri|quad
      |hex|prism|tet|pyramid>
     (required)  element type to be created: line | tri | quad | hex |
     prism | tet | pyramid

   -o <file name of output mesh>,  --mesh-output-file <file name of output
      mesh>
     (required)  the name of the file the mesh will be written to

   --lx <real>
     length of a domain in x direction

   --ly <real>
     length of a domain in y direction

   --lz <real>
     length of a domain in z direction

   --nx <integer>
     the number of subdivision in x direction

   --ny <integer>
     the number of subdivision in y direction

   --nz <integer>
     the number of subdivision in z direction

   --dx0 <real>
     initial cell length in x direction

   --dy0 <real>
     initial cell length in y direction

   --dz0 <real>
     initial cell length in z direction

   --dx-max <real>
     maximum cell length in x direction

   --dy-max <real>
     maximum cell length in y direction

   --dz-max <real>
     maximum cell length in z direction

   --mx <real>
     multiplier in x direction

   --my <real>
     multiplier in y direction

   --mz <real>
     multiplier in z direction

   --ox <real>
     mesh origin (lower left corner) in x direction

   --oy <real>
     mesh origin (lower left corner) in y direction

   --oz <real>
     mesh origin (lower left corner) in z direction


   Structured mesh generator.

   The documentation is available at
   https://docs.opengeosys.org/docs/tools/meshing/structured-mesh-generatio
   n.

   

   OpenGeoSys-6 software, version 6.5.3-407-g9a1da6a18.

   Copyright (c) 2012-2024, OpenGeoSys Community
   (http://www.opengeosys.org)

Quadrilateral mesh

In this step a simple rectangular mesh of size 3 by 4 units, with 4 and 5 cells respectively will be created and written to a quad.vtu file in folder output.
Please follow the folder and file names provided in each step, as some files will be reused later.

In order to select the type of the elements, that the mesh is built with you use the flag -e followed by one of the available mesh element types.
In this section a mesh of quadrilateral elements is created, so the command will be as follows:

generateStructuredMesh -e quad

this however, is not a complete command.
There are more required parameters, next the name of the output file has to be provided.
It is done using the flag -o followed by a file name.
If the mesh file should be written to a folder, it can be provided as path relative to the place from where the command is run.
For example:

generateStructuredMesh -e quad -o folder1/folder2/mesh_file_name.vtu

There are two more groups of parameters that need to be provided: length of the mesh in each direction and number of cells in each direction.
The length parameters are indicated by --l$direction $length, where $direction is the direction indicated by x, y or z and $length is length of the mesh.
In order to create mesh with length of 3 in x-direction and 4 in y-direction.
The number of cells is provided in the same convention with parameter --n: --n$direction $length.
In order to create 4 cells in x-direction and 5 cells in y-direction.
The complete command can be written as follows:

generateStructuredMesh -e quad -o {out_dir}/quads.vtu --lx 3 --ly 4 --nx 4 --ny 5

The created mesh can be viewed with any software supporting .vtu-file format and should look as in the following figure.

Let’s create the mesh with the above command and visualize it:

! generateStructuredMesh -e quad -o {out_dir}/quad.vtu --lx 3 --ly 4 --nx 4 --ny 5… (click to toggle)

! generateStructuredMesh -e quad -o {out_dir}/quad.vtu --lx 3 --ly 4 --nx 4 --ny 5
assert _exit_code == 0
mesh = pv.read(f"{out_dir}/quad.vtu")
show_mesh(mesh)

[2024-12-20 13:28:28.302] [ogs] [[32minfo[m] Mesh created: 30 nodes, 20 elements.

png

Triangle mesh

In this step a mesh using triangular cells will be created.
The size and number of cells will remain the same as in previous example.
The produced mesh file will be placed in output folder in triangle.vtu file.

Compared to the command from previous step, only element type and output file name changes, therefore we can use it:

generateStructuredMesh -e quad -o out/quads.vtu --lx 3 --ly 4 --nx 4 --ny 5

after changing element type quad to tri and file name quad.vtu to triangle.vtu:

! generateStructuredMesh -e tri -o {out_dir}/triangle.vtu --lx 3 --ly 4 --nx 4 --ny 5
assert _exit_code == 0

[2024-12-20 13:28:28.657] [ogs] [[32minfo[m] Mesh created: 30 nodes, 40 elements.

Running this command should create following mesh:

mesh = pv.read(f"{out_dir}/triangle.vtu")
show_mesh(mesh)

png

Shifting the reference point of the mesh

By adding parameter --o$direction $value, the reference point of the created mesh can be shifted.

In this example, the mesh from the step Quadrilateral mesh, will be created again but with the reference point shifted from (0,0) to (2,5) (in (x,y) format).

In order to do this, the final command from that step will be expanded by --ox 2 --oy 5:

! generateStructuredMesh -e quad -o {out_dir}/quads.vtu --lx 3 --ly 4 --nx 4 --ny 5 --ox 2 --oy 5
assert _exit_code == 0

[2024-12-20 13:28:28.986] [ogs] [[32minfo[m] Mesh created: 30 nodes, 20 elements.

Running this command should create following mesh:

mesh = pv.read(f"{out_dir}/quads.vtu")
show_mesh(mesh)

png

Compare the values on X and Y axis with the ones in the figure in step Quadrilateral mesh.

Quadrilateral graded meshes

…with automatically computed cell sizes

In this step a quadrilateral mesh with size of 7 in x-direction and 4 in y-direction (10 and 5 cells respectively) will be created.
However, this time an additional parameter will be passed.
It will change the sizes of the cells from equal to increasingly progressing.

By adding --mx 1.2, each next cell (in positive direction alongside x-axis) will be 1.2 times bigger than the previous one.
Value below zero will make cells smaller than the previous ones.

The command to create this mesh:

! generateStructuredMesh -e quad -o {out_dir}/quads_graded.vtu --lx 7 --ly 4 --nx 10 --ny 5 \… (click to toggle)

! generateStructuredMesh -e quad -o {out_dir}/quads_graded.vtu --lx 7 --ly 4 --nx 10 --ny 5 \
    --mx 1.2
assert _exit_code == 0

[2024-12-20 13:28:29.321] [ogs] [[32minfo[m] Mesh created: 66 nodes, 50 elements.

Note that the name of the output file, values in --l and --n were changed compared to the command in Quadrilateral mesh step.
The \ allows to break the line and split command into more than one line to improve its readability.
Running it produces following mesh:

mesh = pv.read(f"{out_dir}/quads_graded.vtu")
show_mesh(mesh)

png

…with automatically computed cell sizes with given initial cell size

The command from previous step can be expanded with parameter --d$direction0 $value to define initial cell size, alongside $direction (x, y, or z).
Let’s say that the initial cell size is 1.
For this case following parameter should be appended: --dx0 1.
Expanded command from previous step:

! generateStructuredMesh -e quad -o {out_dir}/quads_graded_init_size.vtu --lx 7 --ly 4 --ny 5 \… (click to toggle)

! generateStructuredMesh -e quad -o {out_dir}/quads_graded_init_size.vtu --lx 7 --ly 4 --ny 5 \
    --mx 1.2 --dx0 1
assert _exit_code == 0

[2024-12-20 13:28:29.679] [ogs] [[32minfo[m] Mesh created: 30 nodes, 20 elements.

Note, that output file name has been changed, too.

The mesh produced by that command can be seen in this figure:

mesh = pv.read(f"{out_dir}/quads_graded_init_size.vtu")
show_mesh(mesh)

png

…with automatically computed cell sizes with given initial cell size and maximum cell size

In many cases it may be useful to set a maximal cell size. That can be achieved via the --d$direction-max $value parameter.

The command from the previous step can be expanded by: --dx-max 2:

! generateStructuredMesh -e quad -o {out_dir}/quads_graded_max_size.vtu --lx 7 --ly 4 --ny 5 \… (click to toggle)

! generateStructuredMesh -e quad -o {out_dir}/quads_graded_max_size.vtu --lx 7 --ly 4 --ny 5 \
    --mx 1.2 --dx0 1 --dx-max 2
assert _exit_code == 0

[2024-12-20 13:28:30.027] [ogs] [[32minfo[m] Mesh created: 30 nodes, 20 elements.

Note that the length of mesh and numbers of cells has been modified compared to the previous step in order to improve visibility of the effect of the maximum cell size parameter.

The mesh produced by that command can be seen in this figure:

mesh = pv.read(f"{out_dir}/quads_graded_max_size.vtu")
show_mesh(mesh)

png

Structured mesh in 3D

In previous examples only 2D meshes were created.
createStructuredMesh is capable of creating 3D meshes, too.

In order to create one, the element type has to be changed to a 3D one and parameters --lz and --nz need to be passed additionally.
In 3D, the quadrilateral element type has to be replaced by hexahedral one.

! generateStructuredMesh -e hex -o {out_dir}/hexes.vtu --lx 7 --ly 4 --lz 3 --nx 10 --ny 5 --nz 3
assert _exit_code == 0

[2024-12-20 13:28:30.374] [ogs] [[32minfo[m] Mesh created: 264 nodes, 150 elements.

The command above will create the following mesh:

mesh = pv.read(f"{out_dir}/hexes.vtu")… (click to toggle)

mesh = pv.read(f"{out_dir}/hexes.vtu")
plotter = pv.Plotter()
plotter.add_mesh(mesh, show_edges=True)
plotter.window_size = [800, 400]
plotter.add_axes()
plotter.show()

png

Surface extraction

The following was mainly taken from https://www.opengeosys.org/docs/tools/meshing-submeshes/submeshes/

For Finite Element simulations, we need to specify subdomains e.g. for boundary conditions.
One possibility is to precompute the subdomains as meshes by extracting them from our domain.
This way the subdomains additionally contain information to identify the corresponding bulk mesh entities like nodes, elements, and faces of elements.

We can extract subdomains using the OGS tool ExtractBoundary. Type ! ExtractBoundary --help for a basic documentation.

! ExtractBoundary --help
assert _exit_code == 0

USAGE: 

   ExtractBoundary  [--ascii-output] [-o <file name of output mesh>] -i
                    <file name of input mesh> [--] [--version] [-h]


Where: 

   --ascii-output
     If the switch is set use ascii instead of binary format for data in
     the vtu output.

   -o <file name of output mesh>,  --mesh-output-file <file name of output
      mesh>
     the name of the file the surface mesh should be written to

   -i <file name of input mesh>,  --mesh-input-file <file name of input
      mesh>
     (required)  the name of the file containing the input mesh

   --,  --ignore_rest
     Ignores the rest of the labeled arguments following this flag.

   --version
     Displays version information and exits.

   -h,  --help
     Displays usage information and exits.


   Tool extracts the boundary of the given mesh. The documentation is
   available at
   https://docs.opengeosys.org/docs/tools/meshing-submeshes/extract-boundar
   y.

   

   OpenGeoSys-6 software, version 6.5.3-407-g9a1da6a18.

   Copyright (c) 2012-2024, OpenGeoSys Community
   (http://www.opengeosys.org)

In 2D

We start with extracting the boundary of a 2D mesh.
Let’s remind ourselves of the mesh by visualizing it once more.

mesh = pv.read(f"{out_dir}/quads.vtu")
show_mesh(mesh)

png

From above we can see the required arguments for using the ExtractBoundary tool.
We pass in the mesh to extract the boundary from quads.vtu and save the boundary to quads_boundary.vtu.

! ExtractBoundary -i "{out_dir}/quads.vtu" -o "{out_dir}/quads_boundary.vtu"
assert _exit_code == 0

[2024-12-20 13:28:31.019] [ogs] [[32minfo[m] Mesh read: 30 nodes, 20 elements.
[2024-12-20 13:28:31.020] [ogs] [[32minfo[m] 0 property vectors copied, 0 vectors skipped.
[2024-12-20 13:28:31.020] [ogs] [[32minfo[m] Created surface mesh: 18 nodes, 18 elements.

Visualize the output mesh:

boundary_mesh = pv.read(f"{out_dir}/quads_boundary.vtu")
show_mesh(boundary_mesh)

png

As shown above, the boundary of a 2D mesh is just a set of lines.

In 3D

We can use the same command to extract the boundary of a 3D mesh.

# Using a previously generated mesh… (click to toggle)

# Using a previously generated mesh
! ExtractBoundary -i "{out_dir}/hexes.vtu" -o "{out_dir}/hexes_boundary.vtu"
assert _exit_code == 0

[2024-12-20 13:28:31.376] [ogs] [[32minfo[m] Mesh read: 264 nodes, 150 elements.
[2024-12-20 13:28:31.377] [ogs] [[32minfo[m] 0 property vectors copied, 0 vectors skipped.
[2024-12-20 13:28:31.377] [ogs] [[32minfo[m] Created surface mesh: 192 nodes, 190 elements.

And visualize it:

boundary_mesh = pv.read(f"{out_dir}/hexes_boundary.vtu")
show_mesh_3D(boundary_mesh.shrink(0.8))

png

Mesh clipping

The tool removeMeshElements removes those elements from a given input mesh that fulfill a user specified criterion.
The resulting mesh will be written to the specified output file.
Type ! removeMeshElements --help for a basic documentation.

! removeMeshElements --help
assert _exit_code == 0

USAGE: 

   removeMeshElements  -i <file name of input mesh> -o <file name of output
                       mesh> [--inside] [--outside] [--max-value <number>]
                       [--min-value <number>] [--property-value <number>]
                       ... [-n <string>] [-t <point|line|tri|quad|hex|prism
                       |tet|pyramid>] ... [-z] [--x-min <value>] [--x-max
                       <value>] [--y-min <value>] [--y-max <value>]
                       [--z-min <value>] [--z-max <value>] [--invert] [--]
                       [--version] [-h]


Where: 

   -i <file name of input mesh>,  --mesh-input-file <file name of input
      mesh>
     (required)  the name of the file containing the input mesh

   -o <file name of output mesh>,  --mesh-output-file <file name of output
      mesh>
     (required)  the name of the file the mesh will be written to

   --inside
     remove all elements inside the given property range

   --outside
     remove all elements outside the given property range

   --max-value <number>
     maximum value of range for selected property

   --min-value <number>
     minimum value of range for selected property

   --property-value <number>  (accepted multiple times)
     value of selected property to be removed

   -n <string>,  --property-name <string>
     name of property in the mesh

   -t <point|line|tri|quad|hex|prism|tet|pyramid>,  --element-type <point
      |line|tri|quad|hex|prism|tet|pyramid>  (accepted multiple times)
     element type to be removed: point | line | tri | quad | hex | prism |
     tet | pyramid

   -z,  --zero-volume
     remove zero volume elements

   --x-min <value>
     smallest allowed extent in x-dimension

   --x-max <value>
     largest allowed extent in x-dimension

   --y-min <value>
     smallest allowed extent in y-dimension

   --y-max <value>
     largest allowed extent in y-dimension

   --z-min <value>
     smallest allowed extent in z-dimension

   --z-max <value>
     largest allowed extent in z-dimension

   --invert
     inverts the specified bounding box

   --,  --ignore_rest
     Ignores the rest of the labeled arguments following this flag.

   --version
     Displays version information and exits.

   -h,  --help
     Displays usage information and exits.


   Removes mesh elements based on element type, element volume, scalar
   arrays, or bounding box . The documentation is available at
   https://docs.opengeosys.org/docs/tools/meshing/remove-mesh-elements.

   

   OpenGeoSys-6 software, version 6.5.3-407-g9a1da6a18.

   Copyright (c) 2012-2024, OpenGeoSys Community
   (http://www.opengeosys.org)

The user can choose between 4 different removal criteria:

Remove elements by assigned properties, for instance material ids.
Remove elements by element type, for instance remove line elements.
Remove elements that have zero volume.
Remove elements by axis aligned bounding box criterion.

Let’s load the boundary mesh that we have extracted in the last step:

boundary_mesh = pv.read(f"{out_dir}/hexes_boundary.vtu")… (click to toggle)

boundary_mesh = pv.read(f"{out_dir}/hexes_boundary.vtu")
show_mesh_3D(boundary_mesh.shrink(0.8))

xmin, xmax, ymin, ymax, zmin, zmax = boundary_mesh.bounds
boundary_mesh.bounds

png

(0.0, 7.0, 0.0, 4.0, 0.0, 3.0)

We remove all mesh elements except for the ones at the minimum z-value i.e. the plane at the minimum z-value.

! removeMeshElements -i "{out_dir}/hexes_boundary.vtu" -o "{out_dir}/hexes_zmin.vtu" --z-min {zmin + 1e-8}… (click to toggle)

! removeMeshElements -i "{out_dir}/hexes_boundary.vtu" -o "{out_dir}/hexes_zmin.vtu" --z-min {zmin + 1e-8}
assert _exit_code == 0
mesh = pv.read(f"{out_dir}/hexes_zmin.vtu")
show_mesh_3D(mesh.shrink(0.8))

[2024-12-20 13:28:32.002] [ogs] [[32minfo[m] Mesh read: 192 nodes, 190 elements.
[2024-12-20 13:28:32.002] [ogs] [[32minfo[m] Bounding box of "hexes_boundary" is
x = [0.000000,7.000000]
y = [0.000000,4.000000]
z = [0.000000,3.000000]
[2024-12-20 13:28:32.002] [ogs] [[32minfo[m] 140 elements found.
[2024-12-20 13:28:32.002] [ogs] [[32minfo[m] Removing total 140 elements...
[2024-12-20 13:28:32.002] [ogs] [[32minfo[m] 50 elements remain in mesh.
[2024-12-20 13:28:32.002] [ogs] [[32minfo[m] Removing total 126 nodes...

png

We can also remove all other bounding planes and plot them.

! removeMeshElements -i "{out_dir}/hexes_boundary.vtu" -o "{out_dir}/hexes_zmax.vtu" --z-max {zmax - 1e-8}… (click to toggle)

! removeMeshElements -i "{out_dir}/hexes_boundary.vtu" -o "{out_dir}/hexes_zmax.vtu" --z-max {zmax - 1e-8}
assert _exit_code == 0
! removeMeshElements -i "{out_dir}/hexes_boundary.vtu" -o "{out_dir}/hexes_xmin.vtu" --x-min {xmin + 1e-8}
assert _exit_code == 0
! removeMeshElements -i "{out_dir}/hexes_boundary.vtu" -o "{out_dir}/hexes_xmax.vtu" --x-max {xmax - 1e-8}
assert _exit_code == 0
! removeMeshElements -i "{out_dir}/hexes_boundary.vtu" -o "{out_dir}/hexes_ymin.vtu" --y-min {ymin + 1e-8}
assert _exit_code == 0
! removeMeshElements -i "{out_dir}/hexes_boundary.vtu" -o "{out_dir}/hexes_ymax.vtu" --y-max {ymax - 1e-8}
assert _exit_code == 0

[2024-12-20 13:28:32.318] [ogs] [[32minfo[m] Mesh read: 192 nodes, 190 elements.
[2024-12-20 13:28:32.318] [ogs] [[32minfo[m] Bounding box of "hexes_boundary" is
x = [0.000000,7.000000]
y = [0.000000,4.000000]
z = [0.000000,3.000000]
[2024-12-20 13:28:32.318] [ogs] [[32minfo[m] 140 elements found.
[2024-12-20 13:28:32.318] [ogs] [[32minfo[m] Removing total 140 elements...
[2024-12-20 13:28:32.318] [ogs] [[32minfo[m] 50 elements remain in mesh.
[2024-12-20 13:28:32.318] [ogs] [[32minfo[m] Removing total 126 nodes...

[2024-12-20 13:28:32.517] [ogs] [[32minfo[m] Mesh read: 192 nodes, 190 elements.
[2024-12-20 13:28:32.517] [ogs] [[32minfo[m] Bounding box of "hexes_boundary" is
x = [0.000000,7.000000]
y = [0.000000,4.000000]
z = [0.000000,3.000000]
[2024-12-20 13:28:32.517] [ogs] [[32minfo[m] 175 elements found.
[2024-12-20 13:28:32.517] [ogs] [[32minfo[m] Removing total 175 elements...
[2024-12-20 13:28:32.517] [ogs] [[32minfo[m] 15 elements remain in mesh.
[2024-12-20 13:28:32.517] [ogs] [[32minfo[m] Removing total 168 nodes...

[2024-12-20 13:28:32.719] [ogs] [[32minfo[m] Mesh read: 192 nodes, 190 elements.
[2024-12-20 13:28:32.719] [ogs] [[32minfo[m] Bounding box of "hexes_boundary" is
x = [0.000000,7.000000]
y = [0.000000,4.000000]
z = [0.000000,3.000000]
[2024-12-20 13:28:32.720] [ogs] [[32minfo[m] 175 elements found.
[2024-12-20 13:28:32.720] [ogs] [[32minfo[m] Removing total 175 elements...
[2024-12-20 13:28:32.720] [ogs] [[32minfo[m] 15 elements remain in mesh.
[2024-12-20 13:28:32.720] [ogs] [[32minfo[m] Removing total 168 nodes...

[2024-12-20 13:28:32.922] [ogs] [[32minfo[m] Mesh read: 192 nodes, 190 elements.
[2024-12-20 13:28:32.922] [ogs] [[32minfo[m] Bounding box of "hexes_boundary" is
x = [0.000000,7.000000]
y = [0.000000,4.000000]
z = [0.000000,3.000000]
[2024-12-20 13:28:32.922] [ogs] [[32minfo[m] 160 elements found.
[2024-12-20 13:28:32.922] [ogs] [[32minfo[m] Removing total 160 elements...
[2024-12-20 13:28:32.922] [ogs] [[32minfo[m] 30 elements remain in mesh.
[2024-12-20 13:28:32.922] [ogs] [[32minfo[m] Removing total 148 nodes...

[2024-12-20 13:28:33.122] [ogs] [[32minfo[m] Mesh read: 192 nodes, 190 elements.
[2024-12-20 13:28:33.122] [ogs] [[32minfo[m] Bounding box of "hexes_boundary" is
x = [0.000000,7.000000]
y = [0.000000,4.000000]
z = [0.000000,3.000000]
[2024-12-20 13:28:33.122] [ogs] [[32minfo[m] 160 elements found.
[2024-12-20 13:28:33.122] [ogs] [[32minfo[m] Removing total 160 elements...
[2024-12-20 13:28:33.122] [ogs] [[32minfo[m] 30 elements remain in mesh.
[2024-12-20 13:28:33.122] [ogs] [[32minfo[m] Removing total 148 nodes...

plotter = pv.Plotter()… (click to toggle)

plotter = pv.Plotter()
plotter.add_mesh(pv.read(f"{out_dir}/hexes_xmin.vtu").shrink(0.8), lighting=False, color="red")
plotter.add_mesh(pv.read(f"{out_dir}/hexes_xmax.vtu").shrink(0.8), lighting=False, color="orange")
plotter.add_mesh(pv.read(f"{out_dir}/hexes_ymin.vtu").shrink(0.8), lighting=False, color="green")
plotter.add_mesh(pv.read(f"{out_dir}/hexes_ymax.vtu").shrink(0.8), lighting=False, color="yellow")
plotter.add_mesh(pv.read(f"{out_dir}/hexes_zmin.vtu").shrink(0.8), lighting=False, color="blue")
plotter.add_mesh(pv.read(f"{out_dir}/hexes_zmax.vtu").shrink(0.8), lighting=False, color="cyan")
plotter.add_axes()
plotter.window_size = [800,400]
plotter.show()

png

The special case: generating 0D boundary and source term meshes for OGS

At the moment (Sep 2024) the ExtractBoundary tool unfortunately does not support extracting the zero dimensional boundary of a one dimensional mesh. To create such zero dimensional meshes suitable for source terms and boundary conditions in OGS the following can be done, among others.

First, create a gml file with the following content. Change point coordinates to your needs.

<?xml version="1.0" encoding="ISO-8859-1"?>
<?xml-stylesheet type="text/xsl" href="OpenGeoSysGLI.xsl"?>
<OpenGeoSysGLI xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ogs="http://www.opengeosys.org">
    <name>geometry</name>
    <points>
        <point id="0" x="0.0"  y="0.0" z="0.0" name="source"/>
    </points>
</OpenGeoSysGLI>

Then, run constructMeshesFromGeometry, which will produce a mesh consisting of the single point and set up for use with OGS:

constructMeshesFromGeometry -m SIMULATION_DOMAIN.vtu -g THE_GML_FILE_YOU_JUST_CREATED.gml

Exercises

Exercise 1

Create a 2D mesh consisting of triangles.
There shall be 7 cells in x direction and 6 cells in y direction.
The mesh shall be 3 units wide and 4 units tall.
Its origin shall be at (5, 1).
The result can be visualized with show_mesh().

# Put your solution here

# Visualize your result
mesh_name = "" # Put your mesh name inside the ""
mesh = pv.read(f"out/{mesh_name}")
show_mesh(mesh)

Exercise 2

Consider the mesh from the beginning of the previous step:

! generateStructuredMesh -e quad -o {out_dir}/quads.vtu --lx 3 --ly 4 --nx 4 --ny 5

Remove all elements from this mesh that are not at the mesh boundary.

# Put your solution here

# Visualize your result
mesh_name = "" # Put your mesh name inside the brackets
mesh = pv.read(f"{out_dir}/{mesh_name}")
show_mesh(mesh)

Solutions to the exercises

Exercise 1

! generateStructuredMesh -e tri -o {out_dir}/ex1_triangle.vtu --lx 3 --ly 4 --nx 7 --ny 6 --ox 5 --oy 1… (click to toggle)

! generateStructuredMesh -e tri -o {out_dir}/ex1_triangle.vtu --lx 3 --ly 4 --nx 7 --ny 6 --ox 5 --oy 1
assert _exit_code == 0

# Visualize the result
mesh_name = "ex1_triangle.vtu" # Put your mesh name inside the brackets
mesh = pv.read(f"{out_dir}/{mesh_name}")
show_mesh(mesh)

[2024-12-20 13:28:33.449] [ogs] [[32minfo[m] Mesh created: 56 nodes, 84 elements.

png

Exercise 2

# Extract boundary… (click to toggle)

# Extract boundary
! ExtractBoundary -i "{out_dir}/quads.vtu" -o "{out_dir}/ex2quads_boundary.vtu"
assert _exit_code == 0

# Visualize result
mesh_name = "ex2quads_boundary.vtu" # Put your mesh name inside the brackets
mesh = pv.read(f"{out_dir}/{mesh_name}")
show_mesh(mesh)

[2024-12-20 13:28:33.803] [ogs] [[32minfo[m] Mesh read: 30 nodes, 20 elements.
[2024-12-20 13:28:33.803] [ogs] [[32minfo[m] 0 property vectors copied, 0 vectors skipped.
[2024-12-20 13:28:33.803] [ogs] [[32minfo[m] Created surface mesh: 18 nodes, 18 elements.

png

Which data do the created meshes contain?

Bulk mesh

import pyvista as pv… (click to toggle)

import pyvista as pv

# Created with: ! generateStructuredMesh -e hex -o out/hexes.vtu --lx 7 --ly 4 --lz 3 --nx 10 --ny 5 --nz 3
mesh = pv.read(f"{out_dir}/hexes.vtu")
mesh

UnstructuredGrid	Information
N Cells	150
N Points	264
X Bounds	0.000e+00, 7.000e+00
Y Bounds	0.000e+00, 4.000e+00
Z Bounds	0.000e+00, 3.000e+00
N Arrays	0

From the above cell’s output we see that ‘N Arrays == 0’. This means there is no data contained.

Boundary mesh

mesh = pv.read(f"{out_dir}/hexes_boundary.vtu") # Created by ExtractBoundary tool
mesh

Header

Data Arrays

UnstructuredGrid	Information
N Cells	190
N Points	192
X Bounds	0.000e+00, 7.000e+00
Y Bounds	0.000e+00, 4.000e+00
Z Bounds	0.000e+00, 3.000e+00
N Arrays	3

Name	Field	Type	N Comp	Min	Max
bulk_node_ids	Points	uint64	1	0.000e+00	2.630e+02
bulk_element_ids	Cells	uint64	1	0.000e+00	1.490e+02
bulk_face_ids	Cells	uint64	1	0.000e+00	5.000e+00

mesh = pv.read(f"{out_dir}/hexes_xmax.vtu") # Created by ExtractBoundary tool and processed by removeMeshElements
mesh

Header

Data Arrays

UnstructuredGrid	Information
N Cells	15
N Points	24
X Bounds	7.000e+00, 7.000e+00
Y Bounds	0.000e+00, 4.000e+00
Z Bounds	0.000e+00, 3.000e+00
N Arrays	3

Name	Field	Type	N Comp	Min	Max
bulk_node_ids	Points	uint64	1	1.000e+01	2.630e+02
bulk_element_ids	Cells	uint64	1	9.000e+00	1.490e+02
bulk_face_ids	Cells	uint64	1	2.000e+00	2.000e+00

What is this data used for?

bulk_node_ids identify where the contributions of BCs and STs will go in the global equation system
bulk_element_ids, bulk_face_ids used with some advanced features, e.g. special BCs and surface flux computation

OGS will complain if the needed data is missing for BC or ST meshes.

bulk_mesh = pv.read(f"{out_dir}/hexes.vtu")… (click to toggle)

bulk_mesh = pv.read(f"{out_dir}/hexes.vtu")
bulk_mesh.point_data["bulk_node_ids"] = np.arange(bulk_mesh.n_points)

boundary_mesh = pv.read(f"{out_dir}/hexes_xmax.vtu")

# move the mesh a bit for better visualization
boundary_mesh.points[:,0] += 1

plotter = pv.Plotter()
plotter.add_mesh(bulk_mesh, scalars="bulk_node_ids", show_edges=True, lighting=False)
plotter.add_mesh(boundary_mesh, scalars="bulk_node_ids", show_edges=True, lighting=False)
plotter.window_size = [800,400]
plotter.add_axes()
plotter.show()

png

bulk_mesh = pv.read(f"{out_dir}/hexes.vtu")… (click to toggle)

bulk_mesh = pv.read(f"{out_dir}/hexes.vtu")
bulk_mesh.cell_data["bulk_element_ids"] = np.arange(bulk_mesh.n_cells)

boundary_mesh = pv.read(f"{out_dir}/hexes_xmax.vtu")

# Move the mesh a bit for better visualization
boundary_mesh.points[:,0] += 1

plotter = pv.Plotter()
plotter.add_mesh(bulk_mesh, scalars="bulk_element_ids", show_edges=True, lighting=False)
plotter.add_mesh(boundary_mesh, scalars="bulk_element_ids", show_edges=True, lighting=False)
plotter.window_size = [800,400]
plotter.add_axes()
plotter.show()

png