To continue illustrating this with the example above, you can do the following:
import meshio
import numpy as np
x = np.array([[0,0],[1,0],[0,1],[1,1], [0,2]])
cells = np.array([[0,1,2],[1,2,3], [2,3,4]], dtype=np.int32)
cell_data = np.array([33, 55,11],dtype=np.int32)
cell_data[np.array([1,2],dtype=np.int32)] = 3
mesh_out = meshio.Mesh(points=x, cells={"triangle":cells}, cell_data={"name_to_read":[cell_data]})
meshio.write("test_mesh.xdmf", mesh_out)
from dolfin import *
mesh = Mesh()
with XDMFFile("test_mesh.xdmf") as xdmf:
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
xdmf.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
dx = Measure("dx", domain=mesh, subdomain_data=mf)
print(assemble(1*dx(33)), assemble(1*dx(3)))
as you observe from the output, the mesh function contains each of the subdomains.