[BUG]: creating topology of manifolds makes incorrect assumptions in the dual graph

[nate-sime]

Summarize the issue

Possibly related to #1465.

Creating a mesh of a manifold always assumes that interior (owned) facets will share more than one cell. In the case of manifolds this may not be true as a facet on a process boundary may share more than one cell, whilst also being connected to another cell off process. See, e.g., https://github.com/FEniCS/dolfinx/blob/main/cpp/dolfinx/mesh/graphbuild.cpp#L417-L459.

The attached image shows that the facets (tdim = 0) which share more than one cell (tdim = 1) are not considered for potential "sharing" between neighbouring processes. The pink coloured entities are ghosts.

Thanks for @jorgensd for his help in tracking this one down.

How to reproduce the bug

Consider the following example

Minimal Example (Python)

import dolfinx
import numpy as np
from mpi4py import MPI

def pprint(*msg):
    print(f"[{MPI.COMM_WORLD.rank}]: {', '.join(map(str, msg))}", flush=True)

def print_stats(mesh, name):
    pprint(f"{name} verts {mesh.topology.index_map(0).size_global:,}"
          f" edges {mesh.topology.index_map(1).size_global:,}")

mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 1, 2,
    cell_type=dolfinx.mesh.CellType.quadrilateral)
mesh.topology.create_entities(1)
print_stats(mesh, "mesh")

entities = np.arange(
    mesh.topology.index_map(1).size_local, dtype=np.int32)

mesh_1, _, _ ,_ = dolfinx.mesh.create_submesh(
    mesh, 1, entities)
print_stats(mesh_1, "submesh tdim 1")

with dolfinx.io.XDMFFile(
        mesh_1.comm, "test-mesh.xdmf", "w",
        encoding=dolfinx.io.XDMFFile.Encoding.HDF5) as fi:
    fi.write_mesh(mesh_1)

with dolfinx.io.XDMFFile(
        MPI.COMM_WORLD, "test-mesh.xdmf", "r",
        encoding=dolfinx.io.XDMFFile.Encoding.HDF5) as fi:
    mesh_read = fi.read_mesh()
print_stats(mesh_read, "read mesh")

Output (Python)

$ mpirun -np 1 python3 example.py 
[0]: mesh verts 6 edges 7
[0]: submesh tdim 1 verts 6 edges 7
[0]: read mesh verts 6 edges 7
$ mpirun -np 2 python3 example.py 
[0]: mesh verts 6 edges 7
[1]: mesh verts 6 edges 7
[0]: submesh tdim 1 verts 6 edges 7
[1]: submesh tdim 1 verts 6 edges 7
[0]: read mesh verts 7 edges 7
[1]: read mesh verts 7 edges 7

Version

main branch

DOLFINx git commit

476c61c3b7c743803ea25086ac1dad0ae9b8d7a7

Installation

From source.

Additional information

No response

[jorgensd]

My suggestion would be that we add a GhostMode.manifold option to mesh creation, that is passed to the graph builder. If this mode is used all facets have to be sent to all processors for ownership determination.

this would mean that non-manifold geometries get no penalty, while 2D/3D manifolds will take longer to create (But be correct) with this option turned on. this strategy would be similar to ghostmode.shared_facet, which is almost always only required for DG problems.

If we think that the option should always be on for manifolds (Even if i think that is a bad idea, as not all manifolds has this feature), we could just determine it internally in create_mesh.

[michalhabera]

Possible duplicate of https://github.com/FEniCS/dolfinx/issues/2374?

[garth-wells]

Duplicate of #2374

[chrisrichardson]

Creating a mesh of a manifold always assumes that interior (owned) facets will share more than one cell.

This is a rather confusing statement. In most meshes, we assume that interior facets are attached to exactly two cells. This is also true for non-branching manifold meshes.

In the special case of branching manifold meshes, we have several problems: firstly, the mesh partitioner may struggle, as the dual graph used for partitioning also assumes that an exterior facet is connected to just one cell, and interior facets are connected to two, which is incorrect.

Anyway, if we leave the mesh partitioner aside, it should be possible to create a topology directly with create_topology()

• the final argument boundary_vertices allows you to specify which vertices are on the process boundary - i.e. they could possibly be shared with another process. I think this might allow you to create the topology you want.

[nate-sime]

Thanks for the input. I'll attempt the implementation as you've proposed and perhaps it could be contributed as a demo.

[jorgensd]

• the final argument boundary_vertices allows you to specify which vertices are on the process boundary - i.e. they could possibly be shared with another process. I think this might allow you to create the topology you want

I guess setting all vertices as boundary vertices will have the suggested effect. I still think this could be controlled with a ghostmode in create_mesh

[chrisrichardson]

Although it is not about ghosting. There are various other pathological meshes, e.g. where cells are connected only by a vertex. I think we can support some things, e.g. branching manifolds, in the future, but probably not all cases.