Propagating surface/volume mesh tags from STL/VTM to Firedrake Mesh

[gnikit]

I am wondering if there are any helper functions or already established workflows if you want to

1. Input a geometry represented by a labelled STL.
2. Generate the volumetric mesh with either Netgen or Gmsh
3. Get the tags map that associates regions of the volume mesh with surfaces to programmatically apply boundary conditions.

OR

Load an existing VTM (where the blocks correspond to the boundary regions) directly into Firedrake. The multiblock names/IDs can be fetched with VTK.

Any help would be greatly appreciated.

---

More details

AFAIK both Netgen and Gmsh do not offer a way to map the input labelled STL to the output labelled volume mesh

Netgen

I came close with Netgen fetching the FaceDescriptors() however the 2 surfaces grouped under region5 were split into 2 surfaces again. Also, I doubt that there is a map of the surfnrs and the solid labels in the STL somewhere stored in Netgen.

from netgen.stl import STLGeometry

filename = "cube_labelled.stl"
ngmesh = STLGeometry(filename).GenerateMesh(maxh=0.01)
print(ngmesh.FaceDescriptors())

Outputs

[surfnr = 1, domin = 1, domout = 0, tlosurf = 0, bcprop = 1, bcname = default, domin_sing = 0, domout_sing = 0, colour = (0, 1, 0, 1),
 surfnr = 2, domin = 1, domout = 0, tlosurf = 0, bcprop = 2, bcname = default, domin_sing = 0, domout_sing = 0, colour = (0, 1, 0, 1),
 surfnr = 3, domin = 1, domout = 0, tlosurf = 0, bcprop = 3, bcname = default, domin_sing = 0, domout_sing = 0, colour = (0, 1, 0, 1),
 surfnr = 4, domin = 1, domout = 0, tlosurf = 0, bcprop = 4, bcname = default, domin_sing = 0, domout_sing = 0, colour = (0, 1, 0, 1),
 surfnr = 5, domin = 1, domout = 0, tlosurf = 0, bcprop = 5, bcname = default, domin_sing = 0, domout_sing = 0, colour = (0, 1, 0, 1),
 surfnr = 6, domin = 1, domout = 0, tlosurf = 0, bcprop = 6, bcname = default, domin_sing = 0, domout_sing = 0, colour = (0, 1, 0, 1)]

Gmsh

With Gmsh things are a bit more complicated. One can preserve the input tags using the following script, however, sharp features inside the already classified tags will not be preserved (see figure below)

import gmsh
import numpy as np


gmsh.initialize()
gmsh.merge("cube_labelled.stl")
print("Input STL entities:", gmsh.model.getEntities())

gmsh.model.mesh.createTopology()
surface_tags = [tag for _, tag in gmsh.model.getEntities(2)]
surface_loop = gmsh.model.geo.addSurfaceLoop(surface_tags)
gmsh.model.geo.addVolume([surface_loop])
gmsh.model.geo.synchronize()

# Reclassify surfaces with respect to an angle
# NOTE: this destroys the input file associations
# gmsh.model.mesh.classifySurfaces(np.radians(45), True, True)
gmsh.model.mesh.createGeometry()
for entity in gmsh.model.getEntities():
    if entity[0] < 2:
        continue
    gmsh.model.addPhysicalGroup(entity[0], [entity[1]], entity[1], name=f"Reg_{entity[1]}")
gmsh.model.mesh.generate(3)
gmsh.write("cube_labelled_gmsh.msh")
gmsh.fltk.run()

gmsh.finalize()

Outputs

Enabling gmsh.model.mesh.classifySurfaces(np.radians(45), True, True) will cause the input tag associations to be lost see

Resources

cube_labelled.stl.txt

[ksagiyam]

Firedrake in general does not directly handle mesh files. The Mesh function accepts several types of mesh files as stated in https://www.firedrakeproject.org/firedrake.html#firedrake.mesh.Mesh, but then simply let PETSc create DMPlex objects from those files, and Firedrake works with the DMPlex objects in the rest of the program. So I think the right approach is to have Gmsh or Netgen do the right thing. @UZerbinati might have some suggestions.

[gnikit]

Thanks @ksagiyam for taking the time to reply to this. I think what you are saying is reasonable, but I am not sure how one would be able to retrieve the input geometries region tags in the Netgen/Gmsh/Other mesh.

To try summarise my original question, I was looking for some sort of workflow that

• receives a generic Geometry with regions tagged e.g. inlet, outlet, wall, etc.
• into a Mesher producing a mesh
• usable by Firedrake with the original region tags from the input Geometry somehow available and associated into the Firedrake Mesh.

The end-goal of all this is to programatically be able to set boundary conditions, sources, etc. in a Firedrake model.

Watching again this talk from @UZerbinati I am guessing that this issue might have popped up before since there should not be a direct mapping of an OCC entity's regions to the mesh produced by Netgen (which is consequently fed into Firedrake)