Gmsh import: RuntimeError: Missing entity-vertex connectivity

mesh
1

Dear all,

I am trying to import an gmsh-generated mesh into fenics. However, since fenics can’t handle mixed topologies yet, I need to split the mesh up into different files. I’ve followed this tutorial.

Unfortunately, reading the surface (last line ft = xdmf.read_meshtags(mesh, name="Grid")) raises an Error: RuntimeError: Missing entity-vertex connectivity.

Below is a minimum “working”-example, raising that error. I am not sure if I just messed up the gmsh-generation, so I also attached the *.geo-file. The actual was generated using the GUI.

Thanks in advance.

# meshing.py
import meshio

from dolfinx.io import XDMFFile
from mpi4py import MPI

### Convert from msh to XDMF

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)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

mesh = meshio.read("mesh.msh")

mesh_surf = create_mesh(mesh, "triangle")
meshio.write("mesh_surf.xdmf", mesh_surf)

mesh_vol = create_mesh(mesh, "tetra")
meshio.write("mesh_vol.xdmf", mesh_vol)

### Read XDMF with dolfin

with XDMFFile(MPI.COMM_WORLD, "mesh_vol.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")

mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)

with XDMFFile(MPI.COMM_WORLD, "mesh_surf.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")


//+ mesh.geo
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1, 1, 1};
//+
Physical Surface("s1") = {6};
//+
Physical Surface("s2") = {1};
//+
Physical Surface("s3") = {5};
//+
Physical Surface("s4") = {2};
//+
Physical Surface("s5") = {4};
//+
Physical Surface("s6") = {3};
//+
Physical Volume("dom") = {1};
2

As the error indicates, you are missing the connectivity between facets and vertices.
These can be created by calling:

mesh.topology.create_connectivity(mesh.topology.dim-1, 0)

prior to reading in the meshtag

3

Thank you very much!