I think it should work with HDF5. Try doing this, in serial:
mesh = Mesh("cube.xml")
subdomains = MeshFunction("size_t", mesh, "cube_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "cube_facet_region.xml")
hdf = HDF5File(mesh.mpi_comm(), "file.h5", "w")
hdf.write(mesh, "/mesh")
hdf.write(subdomains, "/subdomains")
hdf.write(boundaries, "/boundaries")
Then in parallel, to read it back in:
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), "file.h5", "r")
hdf.read(mesh, "/mesh", False)
subdomains = CellFunction("size_t", mesh)
hdf.read(subdomains, "/subdomains")
boundaries = FacetFunction("size_t", mesh)
hdf.read(boundaries, "/boundaries")