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

mesh
1

I’m encountering an interesting problem when working (in parallel) on a mesh resulting from local refinement (in serial). Try running the following MWE, first in serial to generate the refined mesh, and then in parallel to trigger a (hard) error:

mpirun -n 1 python3 mwe.py
mpirun -n 2 python3 mwe.py

Any ideas on how to resolve or get around this problem? (I need to work with dS integrals later, hence the shared facet ghost mode.)

The original mesh can be found here https://home.simula.no/~meg/mesh0.h5, while mwe.py reads as:

from dolfin import *

parameters["refinement_algorithm"] = "plaza_with_parent_facets"
parameters["ghost_mode"] = "shared_facet"
#parameters["ghost_mode"] = None # Stops failing with this

info("Reading mesh")
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), "mesh0.h5", "r")
hdf.read(mesh, "/mesh", False)
hdf.close()

V = FunctionSpace(mesh, "CG", 2)
info("Defined V on mesh0")

if MPI.size(MPI.comm_world) == 1:

    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(mesh.mpi_comm(), "mesh1.h5", "w")
    hdf.write(mesh, "/mesh")
    hdf.close()
    
hdf = HDF5File(mesh.mpi_comm(), "mesh1.h5", "r")
hdf.read(mesh, "/mesh", False)
hdf.close()

V = FunctionSpace(mesh, "CG", 2) # This fails
info("Defined V on mesh1")
2 Likes
2

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
3

Thanks Jørgen! Could you say a bit more about the output of the script? In particular,
(1) What does mesh_out show? and
(2) What does “a list of surfaces edges in four facets” mean?

4

mesh_out.xdmf highlight the cells that are connected to “bad_edges”.

Going from all exterior facets (those only connected to one cell, lets call them “surface facets”), we can find all edges that lie on the exterior of the mesh. Let’s call them “surface edges”.
A “surface edge” should only exist in two “surface facets”. If they exist in more than two surface facets, it implies that a set of cells are only connected through an edge or vertex.