Defining meshfunctions for tagged nodes, edges and facets from gmsh mesh (motivated by use in multiphenics)

mesh
1

I would like to create three meshfunctions for nodes, edges and facets of a mesh, from data in a .geo file that was created using gmsh, and converted to .msh. The facets, edges and nodes of the mesh are tagged with integers in the .geo file, and I want the meshfunctions to retain this information. However, so far only the facet tags are retained, and edges and nodes are all tagged with ‘18446744073709551615’, instead of 1,2,3,…

An example .geo file, with tagged subdomains, edges, and nodes:

Mesh.MshFileVersion = 2.0;

lc=5.0e-2;

// Points
Point(1) = {0.0, 0.0, 0.0, lc};
Point(2) = {1.0, 0.0, 0.0, lc};
Point(3) = {1.0, 1.0, 0.0, lc};
Point(4) = {1.0, 2.0, 0.0, lc};
Point(5) = {0.0, 2.0, 0.0, lc};
Point(6) = {0.0, 1.0, 0.0, lc};

// Edges
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,5};
Line(5) = {5,6};
Line(6) = {6,1};

// Interface
Line(7) = {6,3};

// Subdomain boundaries
Curve Loop(8) = {1,2,-7,6};
Curve Loop(9) = {4,5,7,3};

// Subdomains 
Plane Surface(1) = {8};
Plane Surface(2) = {9};

// Physical Points
Physical Point("points",1) = {1,2,3,4,5,6};

// Physical Curves
Physical Curve("lower", 2) = {1,2,6};
Physical Curve("interface", 3) = {7};
Physical Curve("upper", 4) = {3,4,5};

// Physical Surfaces
Physical Surface("lower", 5) = {1};
Physical Surface("upper", 6) = {2};

The code I’m using to output the meshfunctions is:

from dolfin import *
import meshio

msh = meshio.read("mesh.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={'triangle': msh.cells['triangle']}))
meshio.write("nodes.xdmf", meshio.Mesh(points=msh.points, cells={'vertex': msh.cells['vertex']}, cell_data={ \
                             'vertex': {'name_to_read': msh.cell_data['vertex']['gmsh:physical']}}))
meshio.write("edges.xdmf", meshio.Mesh(points=msh.points, cells={'line': msh.cells['line']}, cell_data={ \
                             'line': {'name_to_read': msh.cell_data['line']['gmsh:physical']}}))
meshio.write("facets.xdmf", meshio.Mesh(points=msh.points, cells={'triangle': msh.cells['triangle']}, cell_data={ \
                             'triangle': {'name_to_read': msh.cell_data['triangle']['gmsh:physical']}}))

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

mvc0 = MeshValueCollection("size_t", mesh,0)
with XDMFFile("nodes.xdmf") as infile:
    infile.read(mvc0,'name_to_read')
nodes = cpp.mesh.MeshFunctionSizet(mesh, mvc0)

mvc1 = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("edges.xdmf") as infile:
    infile.read(mvc1, 'name_to_read')
edges = cpp.mesh.MeshFunctionSizet(mesh, mvc1)

mvc2 = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("facets.xdmf") as infile:
    infile.read(mvc2, 'name_to_read')
facets = cpp.mesh.MeshFunctionSizet(mesh, mvc2)

File("nodes.xml") << nodes
File("edges.xml") << edges
File("facets.xml") << facets

I’ve also tried using dolfin-convert .msh to .xml to create the meshfunctions, but this only outputs mesh.xml, mesh-physical-region.xml, and mesh-facet-region.xml; whereas I also want an .xml file containing the node data.

I’m using fenics 2018.1.0, installed on macOS 10.13.6 through conda-forge.

The context to this question is: I’m using gmsh to create a 2D mesh with multiple domains, and using multiphenics to solve different PDEs in different domains, coupled to each other. In doing so, multiphenics defines bool meshfunctions of dimension 0, 1 and 2 in each domain (using MeshRestriction). As I want to use a mesh made in gmsh, as opposed to directly defining the mesh in fenics, the most direct way I can see to construct the bool meshfunctions is through first defining size_t meshfunctions in each subdomain, then converting them to bool. But if anyone has other suggestions or better understanding I’d be grateful!