FEniCS: Read mesh from msh (XDMF) with subdomains

The overall procedure is generating Mesh(), MeshValueCollection objects first, and reading information into them. Then generating MeshFunctionSizet accordingly to define dx, ds... for subdomains.

Step 0. Open corresponding XDMF files; e.g., 

with XDMFFile(MPI.comm_world, "mesh.xdmf") as xdmf_infile:

Step 1. Create Mesh object, using

mesh = dolfin.cpp.mesh.Mesh()

Step 2. Read Mesh from XDMF file, e.g., 

xdmf_infile.read(mesh)

Step 3. Create MeshValueCollection object as container of subdomain information;

mvc_subdomain = dolfin.MeshValueCollection("size_t", mesh, mesh.topology().dim())

Step 4. Read subdomain information,

xdmf_infile.read(mvc_subdomain, "name-to-read")

The "name-to-read" here is the name given in XDMF file generated by meshio.

Step 5. Using step 3, 4 to read boundaries (open another file if necessary):

mvc_boundaries = dolfin.MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
xdmf_infile.read(mvc_boundaries, "name-to-read")

Step 6. Generate MeshFunctionSizet object, this is crucial for the latter solving process:

domains = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc_subdomain)
boundaries = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc_boundaries)

Step 7. Define dx, ds for subdomains and boundaries:

dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

Then we could use dx(1), dx(2)... to define problem in specific sub-domains. The index of sub-domains must be consistent with domains.array(), and same for ds. One other thing notable, the default is to have a unique subdomain with value 0. For example, if domains.array() gives [5,5,5,6,6,6], an error would be raised

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     [email protected]
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to set given (local) rows to identity matrix.
*** Reason:  some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where:   This error was encountered inside PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

To fix the error, just re-set the domain values: domains.set_values([0,0,0,1,1,1]).

Step 8. Give the variational problem for each sub-domain, e.g., for poisson problem:

F = inner(a0 * grad(u), grad(v)) * dx(0) + \
    inner(a1 * grad(u), grad(v)) * dx(1) + \
    inner(g_L, v) * ds(1) - \
    inner(g_R, v) * ds(3) - \
    inner(f, v) * dx(0)  - \
    inner(f, v) * dx(1)