To me this seems like an issue with the original mesh (and therefore the following refinement causes even more issues), as I cannot reproduce it with a built in mesh.
from dolfin import *
parameters["refinement_algorithm"] = "plaza_with_parent_facets"
parameters["ghost_mode"] = "shared_facet"
#parameters["ghost_mode"] = None # Stops failing with this
if MPI.size(MPI.comm_world) == 1:
mesh = UnitCubeMesh(3, 3, 3)
markers = MeshFunction("bool", mesh, 3)
markers[0] = True
#markers.set_all(True) # Stops failing with this
info("Refining mesh")
mesh = refine(mesh, markers, True)
info("Writing mesh")
hdf = HDF5File(MPI.comm_world, "mesh1.h5", "w")
hdf.write(mesh, "/mesh")
hdf.close()
else:
mesh = Mesh()
hdf = HDF5File(MPI.comm_world, "mesh1.h5", "r")
hdf.read(mesh, "/mesh", False)
hdf.close()
V = FunctionSpace(mesh, "CG", 2) # This fails
info("Defined V on mesh1")
@chris and I built a mesh checker to look for topological issues, which is available at:
It checks if the mesh is only connected through vertices or edges, or if an edge exists in more than two surface facets.
Running the script with (needs h5py
, networkx
and meshio
python3 mesh_checker.py --infile mesh0 --topology="mesh/topology" --geometry="mesh/coordinates"
results in a list of surfaces edges in four facets.