Pygmsh Tutorial?

Thank you! Now I can write the two xdmf-files. I still struggle with reading the data back, though. By looking at the link above one should do something like (?):

mesh_file = XDMFFile(self.comm, "mesh.xdmf")
mesh = Mesh()
mesh_file.read(mesh);

mvc = MeshValueCollection("size_t", mesh, 2)
mf_file = XDMFFile(self.comm, "mf.xdmf")
mf_file.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

to read the data back which then could be used as if I had used mshr, that is,

  -snip-

  dx = Measure('dx', domain=mesh, subdomain_data=mf)

 boundary_parts = FacetFunction("size_t", mesh)
 boundary_parts.set_all(0)
 -snip-

However, I get an error:

    mf_file.read(mvc, "name_to_read")
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to find entity in map.
*** Reason:  Error reading MeshValueCollection.
*** Where:   This error was encountered inside HDF5File.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  86617f912d4035520fa988d6af07a19d6069d5ee
*** -------------------------------------------------------------------------