Bug in extracting physical tags from sphere

mesh
1

I want to extract physical tags from my geometry all the way up to zero dimensional entities i.e. points using this script. However, I stumbled upon a strange bug in the process which appears to be present in gmshio.model_to_mesh too.

When it happens:

  1. The geometry is a sphere. Replacing it with a box removes the problem.
  2. Adding physical tags of dimension 1 (i.e. edges). Ordinarily, it should be ignored by gmshio.model_to_mesh

The following MWE produces it:

import gmsh, numpy as np, mpi4py
from dolfinx.io import gmshio

fov_R = 1

gmsh.finalize()
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()

gdim = 3
occ = gmsh.model.occ
fov = occ.addSphere(0, 0, 0, fov_R) # does not work
# fov = occ.addBox(0, 0 ,0, fov_R, fov_R, fov_R) # works
occ.synchronize()

# number all domains
all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

# number all edges
all_edges = gmsh.model.getEntities(gdim - 2)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(gdim)

model_rank = 0
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)

Error message:

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[51], line 35
     32 gmsh.model.mesh.generate(gdim)
     34 model_rank = 0
---> 35 mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)

File ~/miniconda3/envs/fenics-complex-0.9/lib/python3.13/site-packages/dolfinx/io/gmshio.py:245, in model_to_mesh(model, comm, rank, gdim, partitioner, dtype)
    243 # Get mesh geometry and mesh topology for each element
    244 x = extract_geometry(model)
--> 245 topologies = extract_topology_and_markers(model)
    247 # Extract Gmsh cell id, dimension of cell and number of nodes to
    248 # cell for each
    249 num_cell_types = len(topologies.keys())

File ~/miniconda3/envs/fenics-complex-0.9/lib/python3.13/site-packages/dolfinx/io/gmshio.py:139, in extract_topology_and_markers(model, name)
    133 for entity in entities:
    134     # Get cell type, list of cells with given tag and topology
    135     # of tagged cells NOTE: Assumes that each entity only have
    136     # one cell-type,
    137     # i.e. facets of prisms and pyramid meshes are not supported
    138     (entity_types, entity_tags, entity_topologies) = model.mesh.getElements(dim, tag=entity)
--> 139     assert len(entity_types) == 1
    141     # Determine number of local nodes per element to create the
    142     # topology of the elements
    143     properties = model.mesh.getElementProperties(entity_types[0])

AssertionError:
2

There was work on this in Extract edge tags and vertex tags from gmsh and return a NamedTuple by finsberg · Pull Request #3547 · FEniCS/dolfinx · GitHub and Refactor gmshio to use correct naming of subentities by jorgensd · Pull Request #3732 · FEniCS/dolfinx · GitHub

It will be in the next release, but you can easily get a preview of it right now.

A possible suggestion would be to

  1. make a backup of the gmshio.py file in your current installation, just in case you need to go back.
  2. download patches from https://github.com/FEniCS/dolfinx/pull/3732.patch and https://github.com/FEniCS/dolfinx/pull/3547.patch
  3. retain only the parts that change gmshio.py, i.e. ignore changes to demo
  4. apply the simplified patches, in order (3547 first, then the other)
  5. have a look at the changes to the demos in the PR to understand what you should change in your own script.
1 Like
3

Thanks a lot Francesco for your quick reply and detailed instructions to fix it. I willy try and get back if there are issues.

Out of curiosity, is there an expected timeline for the new release?

4

Hopefully during the summer, but we can’t be more precise than that at the moment.

1 Like