Meshio does not read Gmsh physical groups

mesh
18

Hello.
I am also little confused about how to properly retrieve physical names from GMSH mesh to XDMF format for use with FEniCS. I followed the same script as https://fenicsproject.discourse.group/t/transitioning-from-mesh-xml-to-mesh-xdmf-from-dolfin-convert-to-meshio/412/172?u=nareniitr provided earlier for a simple 2D rectangular geometry. I convert .msh file to two xdmf files. When I tried to integrate over area dx and boundary ds, I get correct values. However, when I tried to integrate on just one boundary (e.g. Inlet), the integral value is zero. It appears that subdomain_id is not being recognized. Below is a minimum working example.

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

line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mf2D.xdmf", line_mesh)

triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("mesh2D.xdmf", triangle_mesh)
# fenics imported later to avoid compatibilty issue with h5
from fenics import *
mesh = Mesh()
with XDMFFile("mesh2D.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)

with XDMFFile("mf2D.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

print('num_cells   :', mesh.num_cells())
print('num_vertices:', mesh.num_vertices())
print('cell_type   :', mesh.ufl_cell())

dx=Measure("dx", domain=mesh, subdomain_data=mf)
print(assemble(1*dx))

ds=Measure("ds", domain=mesh, subdomain_data=mf)
print(assemble(1*ds))

ds_1=Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=5)
print(assemble(1*ds_1))

I have also attached geo file I got from GMSH attached here.

//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {1, 0, 0, 1.0};
//+
Point(3) = {1, 1, 0, 1.0};
//+
Point(4) = {0, 1, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {4, 1};
//+
Line(4) = {3, 4};
//+
Curve Loop(1) = {4, 3, 1, 2};
//+
Plane Surface(1) = {1};
//+
Physical Curve("Inlet", 5) = {3};
//+
Physical Curve("Outlet", 6) = {2};
//+
Physical Curve("Top", 7) = {4};
//+
Physical Curve("Bottom", 8) = {1};

//+
MeshSize {1, 2, 3, 4} = 0.05;

Can you please help me understand what is the problem here?