Terribly hacky MVC workaround for pygmsh-meshio-dolfin boundary tagging

meshtagging_mvc.py

from pygmsh.built_in.geometry import Geometry
from pygmsh import generate_mesh
import meshio
import dolfin

geom = Geometry()

lc = .1
p0 = geom.add_point([0, 0, 0], lcar=lc)
p1 = geom.add_point([1, 0, 0], lcar=lc)
p2 = geom.add_point([1, 1, 0], lcar=lc)
p3 = geom.add_point([0, 1, 0], lcar=lc)

l0 = geom.add_line(p0, p1)
l1 = geom.add_line(p1, p2)
l2 = geom.add_line(p2, p3)
l3 = geom.add_line(p3, p0)

ll = geom.add_line_loop(lines=[l0, l1, l2, l3])
ps = geom.add_plane_surface(ll)

# Tag line and surface
geom.add_physical_line(lines=l3, label=444)
geom.add_physical_surface(surfaces=ps, label=456)

print("\n".join(geom._GMSH_CODE))

points, cells, point_data, cell_data, field_data = generate_mesh(geom)

print("Writing 2d mesh for dolfin Mesh")
meshio.write("mesh_2d.xdmf", meshio.Mesh(
    points=points,
    cells={"triangle": cells["triangle"]}))

print("Writing 1d line data for dolfin MVC")
meshio.write("mvc_1d.xdmf", meshio.Mesh(
    points=points,
    cells={"line": cells["line"]},
    cell_data={"line": {"name_to_read": cell_data["line"]["gmsh:physical"]}}
))

mesh_2d = dolfin.Mesh()
with dolfin.XDMFFile("mesh_2d.xdmf") as infile:
    print("Reading 2d mesh into dolfin")
    infile.read(mesh_2d)

mvc = dolfin.MeshValueCollection("size_t", mesh_2d, 1)

with dolfin.XDMFFile("mvc_1d.xdmf") as infile:
    print("Reading 1d line data into dolfin mvc")
    infile.read(mvc, "name_to_read")

print("Constructing MeshFunction from MeshValueCollection")
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh_2d, mvc)