Extracting boundaries of a 2D mesh

There are situations, when you have a 2D domain mesh, but not the corresponding boundary meshes (e.g. when extracting a slice from a 3D model). But you need them to properly set boundary conditions. For those cases ogstools provides a function to generate the individual boundary meshes from the domain mesh or from a continuous boundary mesh.

import ogstools as ot
from ogstools import examples

domain = examples.load_meshseries_THM_2D_PVD()[0]

We can generate the boundary meshes from the given example in the following way and get a dictionary of name and mesh pairs per edge. For details, have a look into the documentation: Meshes. from_mesh().

meshes = ot.Meshes.from_mesh(domain)
for name, mesh in meshes.subdomains.items():
    print(name, mesh)

top UnstructuredGrid (0x7f3d66b3d300)
  N Cells:    83
  N Points:   84
  X Bounds:   -1.400e+03, 8.200e+03
  Y Bounds:   4.064e+01, 8.672e+01
  Z Bounds:   6.700e+03, 6.700e+03
  N Arrays:   0
bottom UnstructuredGrid (0x7f3d66b3eaa0)
  N Cells:    83
  N Points:   84
  X Bounds:   -1.400e+03, 8.200e+03
  Y Bounds:   -1.403e+03, -8.125e+02
  Z Bounds:   6.700e+03, 6.700e+03
  N Arrays:   0
left UnstructuredGrid (0x7f3d66b3fa60)
  N Cells:    44
  N Points:   45
  X Bounds:   -1.400e+03, -1.400e+03
  Y Bounds:   -1.403e+03, 4.186e+01
  Z Bounds:   6.700e+03, 6.700e+03
  N Arrays:   0
right UnstructuredGrid (0x7f3d66b3d3c0)
  N Cells:    44
  N Points:   45
  X Bounds:   8.200e+03, 8.200e+03
  Y Bounds:   -8.125e+02, 6.280e+01
  Z Bounds:   6.700e+03, 6.700e+03
  N Arrays:   0

Let’s display and save them:

fig = meshes.plot()

meshes.save()  # optionally, provide a path

[PosixPath('/tmp/ogstools_root/Meshes/20260323_113703_252682/domain.vtu'), PosixPath('/tmp/ogstools_root/Meshes/20260323_113703_252682/bottom.vtu'), PosixPath('/tmp/ogstools_root/Meshes/20260323_113703_252682/left.vtu'), PosixPath('/tmp/ogstools_root/Meshes/20260323_113703_252682/right.vtu'), PosixPath('/tmp/ogstools_root/Meshes/20260323_113703_252682/top.vtu'), PosixPath('/tmp/ogstools_root/Meshes/20260323_113703_252682/meta.yaml')]

By having the top mesh as a boundary, you can calculate a depth-based water column as an initial conditions or for calculations of integrity criteria

depth = ot.mesh.depth(meshes.domain, meshes["top"])
meshes.domain.point_data["pressure"] = 1000 * 9.81 * depth
fig = ot.plot.contourf(meshes.domain, "pressure")

If you need to model an excavation or similar, the Meshes class provides the useful method remove_material() which removes a specified material from the domain and updates the boundary meshes accordingly. The following example is only for demonstration and is not meant to make practical sense.

x = domain.cell_centers().points[:, 0]
mat_ids = domain["MaterialIDs"]
mat_ids[(mat_ids <= 3) & (x < 0)] = 99
meshes.remove_material(99)
fig = meshes.plot()