Read 2D mesh from XDMFFile

I can't seem to find a way to specify that the input mesh from an XDMF-format should be 2D. Minimal example attached. The mesh that is created with pygmsh (4.3.4) and meshio (2.0.4). However, this example crashes as mesh.geometry()=3.

import meshio
import pygmsh
geometry = pygmsh.built_in.Geometry()
square = geometry.add_rectangle(0, 1, 0, 1, 0)
surface = geometry.add_physical_surface(square.surface,label=1)
mesh_data = pygmsh.generate_mesh(geometry)

points, cells, point_data, cell_data, field_data = mesh_data

meshio.write("mesh.xdmf", meshio.Mesh(points=points,
                                      cells={"triangle": cells["triangle"]}))

from dolfin import *
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

for facet in facets(mesh):
    facet.normal(0)

Writing the file to Legacy xml with the following additional commands and then loading it works.

mesh_data = pygmsh.generate_mesh(geometry, geo_filename="mesh.geo")
import os;
os.system("gmsh -2 mesh.geo")
os.system("dolfin-convert mesh.msh mesh.xml")
mesh = Mesh("mesh.xml")

Comments (5)