FInAT: Dual evaluation for compositions of TensorProductElement + EnrichedElement

[pbrubeck]

Stack

Merge in this order, each PR retargets to main once the one above it lands:

#	PR	what
1	#271	FInAT: select H(div)/H(curl) components with a Delta
2	#268	FInAT: dual evaluate on each sub-element's own points
3	firedrakeproject/firedrake#5289	p-multigrid: remove custom interpolation

Interpolating into the facet space of a hexahedron raised
NotImplementedError: Dual basis not defined for element EnrichedElement,
and matrix-free adjoint interpolation out of a restricted space failed in
gem.unconcatenate. Both come from the same gap: dual_basis is a single weight
tensor on a single point set, and an element whose functionals do not all evaluate
on the same points has no such thing.

An element is a direct sum

Such an element is a direct sum, so give every element the sub_elements it is the
direct sum of — itself, unless it is one — and dual evaluate each on its own points:

• an EnrichedElement is the sum of the elements it enriches;
• a TensorProductElement with a summed factor is the sum of the products of that
  factor's sub-elements with the factors either side of it, since the tensor product
  distributes over the direct sum;
• the H(div)/H(curl), flattened and discontinuous wrappers carry the decomposition of
  whatever they wrap.

dual_basis keeps its exact meaning; it is simply only defined on an element that is
not a sum. A sub-element is an element rather than a bare (Q, x) pair because
contracting it needs its own dual_transformation — a physically mapped
sub-element's basis transformation applies to its own index range, not the enriched
one — and its own get_indices(). dual_evaluation on a sum is then literally
"concatenate my sub-elements' own dual evaluations", and each of those is the
ordinary single-element path.

Why the points are never brought together

The blocks could instead be aligned onto one point set and stacked along a single
point index. They must not be: unconcatenate splits a Concatenate only along an
index that is still free in an assignment, and the point index never is — points are
always contracted by the sum in dual_evaluation. A concatenation over points is a
dead end by construction. It survives into the generated code as a Concatenate
temporary, and the per-block structure that motivated aligning them is destroyed
anyway.

Concatenating over the basis index is exactly what unconcatenate was designed
for; its own docstring works the example through for RTCF7, splitting into blocks
of 8x7 and 7x8 that each keep their own tensor structure. That index stays free in
the assignment, so it splits downstream, and each block is contracted on its own
points with whatever structure it has.

One level at a time

The decomposition goes one level deep and dual_evaluation recurses, so each element
concatenates along the index its own sum occupies. That matters for a tensor
product: its sub-elements own a slice of one factor's basis index and the whole of
every other, which is a sub-box of the basis multiindex rather than a contiguous
range of a flat one. Stacking on that factor's own index is correct wherever the
summed factor sits, whereas flattening straight to the elements that are not sums
would renumber the DoFs unless the summed factor happened to be outermost.

The restricted NCE/NCF tests exercise this several levels deep: the summed factor
sits inside a FlattenedDimensions inside an H(curl) or H(div) wrapper.

Restricted H(div)/H(curl) elements are nodal

EnrichedElement infers nodality from pairwise orthogonality, which only recognises
blocks that map to different components. Restricting an element on a tensor product
cell splits each direction into several blocks that map to the same component, so
the inferred answer was False and restrictions of NCE/NCF were reported as
non-nodal. Restriction selects disjoint subsets of the DoFs, so the restricted
subelements are nodal whenever the original ones are; carry the flag through instead
of re-deriving it, as the tensor product restriction already does.

Testing

check_nodal — applying the dual basis to the primal basis and asserting the
identity — is the mathematical check, and it runs over the enriched, tensor-product
and restricted H(div)/H(curl) cases. The FInAT and gem suites are green (352 passed,
8 skipped), as is
tests/firedrake/multigrid/test_hiptmair.py::test_gmg_hiptmair_hcurl[hexahedron-aij],
which was the failing case that motivated this.

The broader Firedrake suites (multigrid, interpolation) have not been run against
this yet.

[rckirby]

It looks like there are a lot of tests. We're returning a lot of more complicated GEM (that gives fewer flops in the end), so I guess this wants testing in Firedrake as well?

[rckirby]

This seems quite technical, but the tests seem to cover the particular things being added.