Problem converting gmsh to xdmf for fenics

mesh
1

I’m creating a notched domain using pygmsh and, following Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken, I am attempting to convert it to xdmf fo ruse in a fenics code. My code is as follows:

import pygmsh
import meshio
import gmsh
import numpy as np

res = 0.05
geometry = pygmsh.geo.Geometry()
model = geometry.__enter__()

L = 1;
w = 0.5
nl = 0.1
nh = 0.1

points = [model.add_point((-L/2, -w/2, 0), mesh_size=res),
          model.add_point((L/2, -w/2, 0), mesh_size=res),
          model.add_point((L/2, w/2, 0), mesh_size=res),
          model.add_point((nl/2, w/2, 0), mesh_size=res),
          model.add_point((0, w/2-nh, 0), mesh_size=res),
          model.add_point((-nl/2, w/2, 0), mesh_size=res),
          model.add_point((-L/2, w/2, 0), mesh_size=res)]

edges = [model.add_line(points[i], points[i+1])
                 for i in range(-1, len(points)-1)]
notched_bdry = model.add_curve_loop(edges)
interior = model.add_plane_surface(notched_bdry)

model.synchronize()

volume_marker = 6

model.add_physical([notched_bdry], "Interior")
model.add_physical([edges[0]], "Left")
model.add_physical([edges[2]], "Right")
model.add_physical([edges[1],edges[3],edges[4],edges[5],edges[6]], "Neutral")

geometry.generate_mesh(dim=2)

gmsh.write("mesh.msh")
gmsh.clear()
geometry.__exit__()

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

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

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

and this produces the error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [1], line 55
     52 line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
     53 meshio.write("facet_mesh.xdmf", line_mesh)
---> 55 triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
     56 meshio.write("mesh.xdmf", triangle_mesh)

Cell In [1], line 47, in create_mesh(mesh, cell_type, prune_z)
     45 def create_mesh(mesh, cell_type, prune_z=False):
     46     cells = mesh.get_cells_type(cell_type)
---> 47     cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
     48     points = mesh.points[:,:2] if prune_z else mesh.points
     49     out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})

File ~/miniconda3/lib/python3.9/site-packages/meshio/_mesh.py:249, in Mesh.get_cell_data(self, name, cell_type)
    248 def get_cell_data(self, name: str, cell_type: str):
--> 249     return np.concatenate(
    250         [d for c, d in zip(self.cells, self.cell_data[name]) if c.type == cell_type]
    251     )

File <__array_function__ internals>:180, in concatenate(*args, **kwargs)

ValueError: need at least one array to concatenate
2

Simply by inspecting mesh_from_file

In [1]: mesh_from_file
Out[1]: 
<meshio mesh object>
  Number of points: 64
  Number of cells:
    line: 10
    line: 20
    line: 10
    line: 9
    line: 3
    line: 3
    line: 9
  Cell sets: Interior, Left, Right, Neutral, gmsh:bounding_entities
  Point data: gmsh:dim_tags
  Cell data: gmsh:physical, gmsh:geometrical
  Field data: Interior, Left, Right, Neutral

you see that there are no triangular cells in your mesh.
This is because you have not used any Physical Surface in pygmsh, i.e.
adding:

model.add_physical([interior], "Interior Surface")

resolves the issue.