Why are (singular) 2D cells appearing in my 3D mesh?

mesh
1

Hi,

I am using gmsh to generate a 3D domain. The output from gmsh states that there are “no ill-shaped tets in the mesh.” However, when I import the mesh using df.io.gmshio.model_to_mesh, I find that the resulting dolfinx mesh has singular tetrahedra (whose 4 points lie in a 2D plane).

The behavior is somewhat non-deterministic. I find singular cells when I run mpiexec -n 2 python test_singular_cells.py and the same with 4 MPI processes, but I do not find singular cells when I run with 3 MPI processes.

Am I generating this mesh incorrectly, or am I perhaps converting it to a dolfinx mesh incorrectly? I don’t know what relevance it may or may not have, but note that I am having gmsh generate a mesh with boundary faces that match periodically. (See the use of “setPeriodic” in the code below.)

Thank you for any help!

# test_singular_cells.py
from mpi4py import MPI
import numpy as np
import gmsh
import dolfinx as df
########################################################
mpi_comm = MPI.COMM_WORLD
root_rank = 0
this_rank = mpi_comm.rank
########################################################
gdim = 3
dLengths = [1,1,1]
########################################################
tgt_cell_size = .1
if this_rank == root_rank:
    gmsh.initialize()
    gmsh.model.add(f'rectangular_{gdim}D_mesh')
    # Add geometry
    Lx, Ly, Lz = dLengths
    p1 = gmsh.model.occ.addPoint(0, 0, 0)
    p2 = gmsh.model.occ.addPoint(Lx, 0, 0)
    p3 = gmsh.model.occ.addPoint(0, Ly, 0)
    p4 = gmsh.model.occ.addPoint(0, 0, Lz)
    p5 = gmsh.model.occ.addPoint(Lx, Ly, 0)
    p6 = gmsh.model.occ.addPoint(Lx, 0, Lz)
    p7 = gmsh.model.occ.addPoint(0, Ly, Lz)
    p8 = gmsh.model.occ.addPoint(Lx, Ly, Lz)
    lx1 = gmsh.model.occ.addLine(p1, p2)
    lx2 = gmsh.model.occ.addLine(p3, p5)
    lx3 = gmsh.model.occ.addLine(p4, p6)
    lx4 = gmsh.model.occ.addLine(p7, p8)
    ly1 = gmsh.model.occ.addLine(p1,p3)
    ly2 = gmsh.model.occ.addLine(p2,p5)
    ly3 = gmsh.model.occ.addLine(p4,p7)
    ly4 = gmsh.model.occ.addLine(p6,p8)
    lz1 = gmsh.model.occ.addLine(p1,p4)
    lz2 = gmsh.model.occ.addLine(p2,p6)
    lz3 = gmsh.model.occ.addLine(p3,p7)
    lz4 = gmsh.model.occ.addLine(p5,p8)
    loop_xy1 = gmsh.model.occ.addCurveLoop([lx1, ly2, -lx2, -ly1]) # Surface orientation:  inward
    loop_xy2 = gmsh.model.occ.addCurveLoop([lx3, ly4, -lx4, -ly3]) # Surface orientation: outward
    loop_yz1 = gmsh.model.occ.addCurveLoop([ly1, lz3, -ly3, -lz1]) # Surface orientation:  inward
    loop_yz2 = gmsh.model.occ.addCurveLoop([ly2, lz4, -ly4, -lz2]) # Surface orientation: outward
    loop_zx1 = gmsh.model.occ.addCurveLoop([lz1, lx3, -lz2, -lx1]) # Surface orientation:  inward
    loop_zx2 = gmsh.model.occ.addCurveLoop([lz3, lx4, -lz4, -lx2]) # Surface orientation: outward
    surf_xy1 = gmsh.model.occ.addPlaneSurface([loop_xy1])
    surf_xy2 = gmsh.model.occ.addPlaneSurface([loop_xy2])
    surf_yz1 = gmsh.model.occ.addPlaneSurface([loop_yz1])
    surf_yz2 = gmsh.model.occ.addPlaneSurface([loop_yz2])
    surf_zx1 = gmsh.model.occ.addPlaneSurface([loop_zx1])
    surf_zx2 = gmsh.model.occ.addPlaneSurface([loop_zx2])
    sloop = gmsh.model.occ.addSurfaceLoop([surf_xy1, surf_yz1, surf_zx1, surf_yz2, surf_zx2, surf_xy2])
    gmsh.model.occ.addVolume([sloop])
    # Synchronize, before gmsh.model.mesh is used
    gmsh.model.occ.synchronize()
    # Periodicize
    gmsh.model.mesh.setPeriodic(2, [surf_yz2], [surf_yz1], [1,0,0,Lx, 0,1,0,0,  0,0,1,0,  0,0,0,1]) # Affine transformation: A*(vx;vy;vz;1)
    gmsh.model.mesh.setPeriodic(2, [surf_zx2], [surf_zx1], [1,0,0,0, 0,1,0,Ly,  0,0,1,0,  0,0,0,1]) # Affine transformation: A*(vx;vy;vz;1)
    gmsh.model.mesh.setPeriodic(2, [surf_xy2], [surf_xy1], [1,0,0,0,  0,1,0,0, 0,0,1,Lz,  0,0,0,1]) # Affine transformation: A*(vx;vy;vz;1)
    # Now generate the mesh for the gmsh model
    cells = gmsh.model.getEntities(gdim)
    tagPhys = 1+np.arange(len(cells))                                            #gmsh starts count at 1 (not 0)
    gmsh.model.addPhysicalGroup(gdim,tagPhys,name=f'Rectangular {gdim}D Domain') #needed by df.io.gmshio.model_to_mesh
    meshsize_field = gmsh.model.mesh.field.add('MathEval')
    gmsh.model.mesh.field.setString(meshsize_field, 'F', str(tgt_cell_size))
    gmsh.model.mesh.field.setAsBackgroundMesh(meshsize_field)
    gmsh.model.mesh.generate(gdim)
#output: df_mesh, cell_tags, facet_tags
df_mesh, _, _ = df.io.gmshio.model_to_mesh(gmsh.model, mpi_comm, root_rank, gdim, \
                    partitioner = df.mesh.create_cell_partitioner(df.mesh.GhostMode.shared_facet)) 
if this_rank == root_rank:
    gmsh.finalize()
assert(gdim==df_mesh.topology.dim)
########################################################
for idim in range(gdim+1):
    df_mesh.topology.create_entities(idim)
nV = df_mesh.topology.index_map(0).size_local
nC = df_mesh.topology.index_map(gdim).size_local
########################################################
xcoor = df_mesh.geometry.x
df_mesh.topology.create_connectivity(gdim, 0)
for iC in range(nC):
    v_of_c = df_mesh.topology.connectivity(gdim, 0).links(iC)
    if np.linalg.matrix_rank(xcoor[v_of_c])<gdim:
        print(f"PROCESS #{this_rank}: SINGULAR CELL: Vertices {v_of_c} of {nV} local vertices. Coordinates: {xcoor[v_of_c]}.")
2

Please note that the map

a not correct, as the mesh topology and geometry has separate dofmaps.

This is because a topology is always linear (consisting of only vertices), while a mesh geometry can have more nodes than vertices (to support higher order cells with curved boundaries).
You should use:

for iC in range(nC):
    v_of_c = df_mesh.topology.connectivity(gdim, 0).links(iC)
    mesh_nodes = df.cpp.mesh.entities_to_geometry(df_mesh._cpp_object, 0, [v_of_c], False).reshape(-1)
3

Thank you, this was of great help!