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 contiuous 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 (0x7f86fcd561a0)
  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 (0x7f86fcd55c60)
  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 (0x7f86fcd56380)
  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 (0x7f86fcd55b40)
  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/tmpz0lvzwnemeshes/domain.vtu'), PosixPath('/tmp/tmpz0lvzwnemeshes/bottom.vtu'), PosixPath('/tmp/tmpz0lvzwnemeshes/left.vtu'), PosixPath('/tmp/tmpz0lvzwnemeshes/right.vtu'), PosixPath('/tmp/tmpz0lvzwnemeshes/top.vtu')]

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