Defining non-CG1 FunctionSpace in parallel on mesh resulting from local mesh refinement

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.

3 Likes