Fix matrix-free adjoint interpolation into Restricted/EnrichedElement

[pbrubeck]

Stack

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

#	PR	what
1	firedrakeproject/fiat#269	gem: sum factorise independent contractions
2	firedrakeproject/fiat#271	FInAT: select H(div)/H(curl) components with a Delta
3	firedrakeproject/fiat#268	FInAT: dual evaluate on each sub-element's own points
4	#5295	tsfc: dual evaluation against a Cofunction
5	#5289	p-multigrid: remove custom interpolation
6	#5297	tsfc: contract each block after unconcatenating

2 and 6 are a pair: neither is worth anything alone, though 2 is harmless on its
own — it changes how a component is selected, not what is selected.
5 and 6 are siblings on 4, in either order.

Only 4 carries a DROP BEFORE MERGE commit, pinning FIAT to
pbrubeck/fix/dual-enriched, which is now the top of the FIAT stack; 5 and 6
inherit it. Drop it once the FIAT side has landed.

Assembling a two-form Interpolate against a Cofunction — the adjoint of a
matrix-free interpolation matrix, i.e. MatMultTranspose — failed for any target
element with a block-structured dual basis, such as a RestrictedElement:

File "gem/unconcatenate.py", line 103, in find_group
    assert i in free_indices
AssertionError

compile_expression_dual_evaluation contracts the dual evaluation against the dual
coefficient with an IndexSum over the basis indices when the dual argument is a
Cofunction. When the target element's dual evaluation is a gem.Concatenate, that
contraction hides the concatenation index from the subsequent unconcatenate call,
which can only split along indices that are still free in the assignment.

Resolve the concatenation of the dual basis before the basis indices are summed
over, following the same pattern already used for coefficient evaluation in
tsfc.fem. The two unconcatenate calls share an index splitting cache so a given
Concatenate index is always split the same way.

Interpolation points come from the dual point set

firedrakeproject/fiat#268 makes an element whose functionals do not all evaluate on
the same points a direct sum, with a dual basis on each of its sub-elements rather
than one weight tensor on one point set. Reconstructing such a space as a runtime
Quadrature space only ever wanted the points, so it asks for dual_point_set
directly instead of unpacking dual_basis.

Depends on firedrakeproject/fiat#268, which supplies the corresponding gem fix:
find_group must intersect against the assignment variables' free indices, not the
expressions'. Without it the extra unconcatenate call here can recurse on a pair it
cannot split.

Testing

Matrix-free mult/multTranspose were compared against the assembled matrix over
quadrilateral, triangle, hexahedron and prism meshes, CG/RTCE/RTCF, full/facet/interior
restrictions, spectral and fdm variants: 47 cases, all agreeing to machine precision.

The p-multigrid and FDM PR is stacked on top of this one.

🤖 Generated with Claude Code