[achanbour]

Compute and cache properties derived from cell-facet adjacency data on meshes:

• cell_facet_neighbours: maps each cell to its neighbouring cells across each of its local facets
• cell_facet_coord_transforms: stores affine reference-coordinate transforms between neighbouring cells across interior facets
• cell_facet_exterior_mask: marks cell-local facets lying on the exterior domain boundary (distinguishes exterior facets from rank partition boundaries)

[achanbour]

Moved helper methods that operate on the reference element to FIAT, see firedrakeproject/fiat#303

[pbrubeck]

It's hard to see the ultimate goal without a consumer of the attributes added by the PR. Why does the low-level reference-facet data need to be attached to every facet?

[pbrubeck]

The data can be derived from the DMPlex, but DMPlex does not return the arrays in Firedrake’s local facet numbering.

• Neighbouring cells: DMPlex’s facet support gives the incident cells. Firedrake already uses that data in local_facet_number (firedrake/cython/dmcommon.pyx:357)
  and _support_dat (firedrake/mesh.py:1635). To build Cache cell-facet topological data on meshes #5298’s per-cell, per-local-facet arrays, the code still has to map DMPlex point IDs to Firedrake cell numbers and
  match each facet to its FIAT local facet number. The current cell-closure data supplies that ordering; DMPlex support order should not be treated as local facet order.

• Exterior-facet mask: Read the exterior_facets DMLabel. Firedrake creates it with DMPlexMarkBoundaryFaces and already uses it to classify cell-local facets in
  cell_facet_labeling (firedrake/cython/dmcommon.pyx:1479). Don’t infer physical exterior facets from support size alone, since a partition boundary can also have
  incomplete local support.

• Coordinate transforms: DMPlex stores incidence and orientation, not the affine reference-coordinate matrices. Those can be derived from DMPlex’s oriented cones or
  closures plus the reference-cell facet maps, then cached.

[pbrubeck]

The coordinate transforms can be obtained from UFL: CellFacetJacobian (A) CellFacetJacobianInverse (pinv(A)) and CellFacetOrigin (b).

You can project them into

BQ = TensorFunctionSpace("Boundary Quadrature", 0, shape=shape).broken_space()

[achanbour]

@pbrubeck For neighbouring cells, my updated code already uses the DMPlex information. Same for exterior facet masks .

For the cell-facet coordinate transforms, I need to go through the FIAT local numbering (that is store the arrays in FIAT local facet ordering) because my particle tracing algorithms are implemented in the reference cell. The mesh entity collision algorithm figures out which facet was crossed based on the barycentric coordinates of the point at the crossing:

crossed_edge = int(np.argmin(abs(bary_cross[i])))

Which is why the arrays used by the cell-crossing algorithm, namely to figure out which cell is the neighbour across the crossed facet and to map the coordinates of the point to that cell, have to be filled using the same numbering as the crossed facet index (which is FIAT local ordering).

[pbrubeck]

From my understanding you are converting from reference coordinates on the '-' side of a facet to the reference coordinates on the '+' side. Is that correct? And you have to go through the horrible orientation data structures from FIAT.

Firedrake uses these orientation arrays for assembly already, so my suggestion to use UFL was to avoid two places in the code that read them.

However, is not easier to have neighbors communicate to each other in physical coordinates? Say that you compose the transformations: X_ref^- -> X_phys -> X_ref^+.

[achanbour]

Is that correct?

That's right.

Firedrake uses these orientation arrays for assembly already

I'm not entirely clear on this approach. Firstly, in which reference frame do these UFL expressions apply? Are they describing transformations between facet-local coordinates or cell coordinates of points on facets? What I'm computing in my approach is the latter via the chain X_{c1} -> X_{f1} -> X_{f2} -> X_{c2} where f1 and f2 correspond to the same mesh facet.