Unable to optimize higher order mesh

The gmshio.model_to_mesh crashes when it operates on a higher order mesh which was later optimized as well. It works fine if the optimization is done on first order mesh before conversion to higher order. I get the following error (dolfinx 0.9, conda):

Invalid rank, error stack:
internal_Issend(117): MPI_Issend(buf=0x5b02d65a2781, count=1, MPI_BYTE, 3, 1, comm=0xc4000004, request=0x5b02d659a014) failed
internal_Issend(78).: Invalid rank has value 3 but must be nonnegative and less than 1
Abort(738854406) on node 0 (rank 0 in comm 304): application called MPI_Abort(comm=0xC4000004, 738854406) - process 0

Following MWE reproduces the error:

import dolfinx, gmsh, mpi4py
from dolfinx.io import gmshio

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

occ = gmsh.model.occ
gdim = 3

mag_R = 1
mag_L = 1
curv_R = 2

t1 = occ.addTorus(curv_R, 0, mag_L/2, curv_R, mag_R, zAxis=(0, 0, 1))

occ.synchronize()

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

gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 10)
# gmsh.option.setNumber("Mesh.ElementOrder", 2)
gmsh.model.mesh.generate(gdim)
# gmsh.model.mesh.optimize("Netgen") # works
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen") # crashes

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

I inspected the mesh with:

gmsh.write("mwe_gmsh.msh")
import meshio

msh = meshio.read("mwe_gmsh.msh")
print(msh)

and that yields:

<meshio mesh object>
  Number of points: 1028
  Number of cells:
    triangle6: 502
    tetra: 743
  Cell sets: gmsh:bounding_entities
  Point data: gmsh:dim_tags
  Cell data: gmsh:physical, gmsh:geometrical

while removing the optimize yields

<meshio mesh object>
  Number of points: 1560
  Number of cells:
    triangle6: 502
    tetra10: 759
  Cell sets: gmsh:bounding_entities
  Point data: gmsh:dim_tags
  Cell data: gmsh:physical, gmsh:geometrical

meaning that optimize removes the higher order tetrahedra. Please add in

gmsh.model.mesh.setOrder(2)

post calling optimize to make it work. If you would like to pursue this further, I would suggest raising an issue in GMSH: Issues · gmsh / gmsh · GitLab

1 Like