Global dual graph computation for branching meshes

[schnellerhase]

Adds branching mesh support to compute_nonlocal_dual_graph algorithm, i.e. support for arbitrary remote matching.

Main change is the altered data exchange. Before, for every unmatched cell the post office deduced if there was a match or none. This info was encoded in an array with either the cell index of the remote match or $-1$ if unmatched for example

$$ [4, 7, -1, -1, 9]. $$

This data was then communicated back to the sender, which could construct from this information the parallel aware dual graph (local dual graph + ghost edges).

This data protocol has the matching condition encoded into it and can not be extended to the case of arbitrary matching conditions. In particular the necessity to communicate back multiple matches make its not well suited.

The data that needs to be exchanged is the facet to cell connectivity adjacency list (restricted to cross boundary parts). For every facet send to the PO we get back a (possibly empty) list of matched cells. An adjacency list is best described as a combination of size information (offset or link count) and the (continuous) data. This motivates the following altered protocol.

1. Communicate (from postoffice to sender) number of remote matches for every send facet. For the previous example this would be

$$ [1, 1, 0, 0, 1] $$

2. follow up with sending unrolled matched cells data:

$$ [4, 7, 9] $$

This extends easily to multiple matches, for example if facet 1 was not only matched to cell $7$ but also cell $2$, the arrays become

$$ [1, 2, 0, 0, 1] \quad \text{and} \quad [4, 7, 2, 9] $$

Note: Communication steps 1. and 2. are not combinable. Step 1 is necessary to deduce size information for step 2. Therefore one additional communication step is necessary (on neighbourhood level).

Example:

One use case for branching meshes are skeleton meshes. Meshes which are constructed from the lower dimensional entities of another mesh. For example for a regular quadrilateral mesh, taking only the edges and construction a mesh with them we get (cell information):

And in parallel:

process 0	process 1

Fixes #3733.

To be addressed in future PR

• compute_entities_by_key_matching/get_local_indexing need to be updated to for correct connectivity computation for cell->edge and edge->vertex connectivities.

[chrisrichardson]

I wonder if the parameter max_facet_to_cell_links is really useful? As far as I can see, there are only two sensible values: 2 or infinity.

[jorgensd]

I wonder if the parameter max_facet_to_cell_links is really useful? As far as I can see, there are only two sensible values: 2 or infinity.

You can analyze this for certain grids. Many grids might have at max 3-4 connections, which would mean that one can in a quicker way reduce the amount of potential communication in the setup phase.

[chrisrichardson]

How does it reduce the amount of communication in the setup phase?

[schnellerhase]

One such example would be the mesh in the PR description - where max_facet_to_cell_links=4 will allow to rule out locally interior facets prior to communicating them. They are never transferred to the post office for matching.

I agree that this is is not as general as one would like to have it, as it requires some very detailed a-priori knowledge about the mesh. However, we can extend this freely in the future to a facet wise input from the user.

[jorgensd]

One such example would be the mesh in the PR description - where max_facet_to_cell_links=4 will allow to rule out locally interior facets prior to communicating them. They are never transferred to the post office for matching.

I agree that this is is not as general as one would like to have it, as it requires some very detailed a-priori knowledge about the mesh. However, we can extend this freely in the future to a facet wise input from the user.

For 1D grids/graphs one can easily get this information through software such as networkx, and then the input will be very useful indeed

[chrisrichardson]

This line looks wrong...

dolfinx/cpp/dolfinx/mesh/graphbuild.cpp

Line 138 in 1a13864

max_vertices_per_facet = std::int64_t(local_max_vertices_per_facet);

[schnellerhase]

This line looks wrong...

dolfinx/cpp/dolfinx/mesh/graphbuild.cpp

Line 138 in 1a13864

max_vertices_per_facet = std::int64_t(local_max_vertices_per_facet);

Good catch, should be static_cast<std::int32_t>(local_max_vertices_per_facet), right?

[chrisrichardson]

I'm not sure, it is all rather messed up. Probably it can just go as a cast to int64 in the list assignment in line 152 and delete line 138.