Cannot tag spherical surface in 3d mesh

General
1

Hello,
I would like to create a mesh with given by a cube with a spherical hole in it. I can do this with the following code

import meshio
import gmsh
import pygmsh
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("resolution")
args = parser.parse_args()

resolution = (float)(args.resolution)

geometry = pygmsh.occ.Geometry()
model = geometry.__enter__()

L1 = 1
L2 = 1
L3 = 1
r = 0.25
c_r = [0.5, 0.5, 0.5]


box = model.add_box([0, 0, 0], [L1, L2, L3], mesh_size=resolution)
ball = model.add_ball(c_r, r,  mesh_size=resolution)

difference = model.boolean_difference(box, ball)

model.synchronize()
model.add_physical( difference, "difference" )

geometry.generate_mesh(dim=3)
gmsh.write("solution/mesh.msh")
model.__exit__()

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    return out_mesh

mesh_from_file = meshio.read("solution/mesh.msh")

tetrahedron_mesh = create_mesh(mesh_from_file, "tetra", True)
meshio.write("solution/tetrahedron_mesh.xdmf", tetrahedron_mesh)

which works

$clear; clear; rm -r solution; mkdir solution; python3 generate_3dmesh_box_ball.py 0.1

and mesh.msh looks like this:

Screenshot 2025-01-08 at 14.48.14
Screenshot 2025-01-08 at 14.48.141610×1382 393 KB

I can tag the volume given by the difference between the cube and the ball with model.add_physical( difference, "difference" ) which is tagged with name_to_read = 1.

I would like to tag the ball surface, only the surface, not its volume, but I don’t know how to do it. Do you know how to do this ?

Note that I tried multiple solutions: this shows how to tag a sphere with its volume, which is not what I want. This allowed me to extract the surfaces with surfaces = gmsh.model.occ.getEntities(dim=2) but then it is unclear how to tag them: The webpage uses gmsh.model.addPhysicalGroup while in the code above I use model.add_physical, which does not work with argument surfaces[0].

Thank you

2

I would assume there is an easier way, but I’ve been looping over all ‘boundaries’ as suggested in: Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial

So, you’d get something like:

sphere_tag = 99 # Can be used in your FEniCS code later as a facet marker
dim_facet = 2 # for facets in 3D
sphere_boundaries = []
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
    center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
    if np.sum(center_of_mass**2) < radius**2 + 1e-7:
        sphere_boundaries.append(boundary[1])
gmsh.model.addPhysicalGroup(dim_facet, sphere_boundaries, tag=sphere_tag)
4

The problem/confusion here is that the first example refers to pygmsh, which is an interface on top of the gmsh api.

The second example uses the gmsh python api directly.

To give some context regarding the two different ways of using gmsh:

  1. At its inception Pygmsh was a nice way of creating gmsh meshes through Python. It would generate .geo files.
  2. After some time, the developer of the library decided to switch to the gmsh Python api as backend of pygmsh, meaning that pygmsh became a way of calling specific gmsh calls.

Personally, this mean that I got more freedom and a cleaner interface interfacing directly with gmsh, which is the way I would recommend using the software these days.

Long story short; I would stop using pygmsh and use the gmsh Python api directly.

2 Likes
5

I must have made a mistake when trying to apply this solution: by going through it again, it works. I included the lines that you posted in mine, then exported the mesh of triangles with

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=False)
meshio.write("solution/triangle_mesh.xdmf", triangle_mesh)

and when I read back triangle_mesh.xdmf I find all the surfaces as expected :slight_smile:
The ids of these surfaces are assigned in a quite odd way, but one can get them easily with Paraview by opening triangle_mesh.xdmf and checking the color for name_to_read assigned to the different surfaces.

Thank you for your solution!

6

Thank you for the clarification