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
*** -------------------------------------------------------------------------