Select H(div)/H(curl) components with a Delta

[pbrubeck]

Stack

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

#	PR	what
1	#269	gem: sum factorise independent contractions
2	#268	FInAT: EnrichedElement dual basis
3	#270	FInAT: select H(div)/H(curl) components with a Delta
4	firedrakeproject/firedrake#5295	tsfc: dual evaluation against a Cofunction
5	firedrakeproject/firedrake#5289	p-multigrid: remove custom interpolation
6	firedrakeproject/firedrake#5297	tsfc: contract each block after unconcatenating

3 and 6 are a pair and must land together: neither does anything alone.
5 and 6 are siblings on 4, in either order.

The Firedrake PRs carry a DROP BEFORE MERGE commit pinning FIAT to the
paired branch; drop it once the FIAT side has landed.

An H(div) or H(curl) element is a sum of blocks, each mapping the value of a
wrapped tensor product element into one component of a vector. The transformers
expressed that by padding a ListTensor with zeros:

return lambda v: [gem.Zero(), gem.Zero(), v]

Contracting two such blocks over the value index cannot cancel the pairs that map
to different components, because the zeros are entries of a tensor rather than a
property of the expression. So an interpolation between two H(div) spaces keeps
every pair of blocks in the kernel, including the ones that are identically zero.

Select the component with a Delta instead:

return lambda v, zeta: gem.Product(gem.Delta(zeta, 2), v)

A contraction of two blocks then carries a product of Deltas, which delta
elimination cancels unless the blocks map to the same component. The transform now
takes the value index and returns the component it selects, rather than returning a
list of components to be indexed.

Paired with firedrakeproject/firedrake#5297, which contracts each block once
unconcatenate has separated them — that is where the Deltas meet pairwise.
Neither PR does anything on its own.

Effect

Coarse to fine interpolation, as p-multigrid prolongation. Kernel stack
temporaries and wall clock:

	temporaries	time
NCF 1→4	832K → 68K	3.07×
NCE 1→4	307K → 110K	2.23×
NCE 1→3	131K → 61K	2.23×
NCF 1→3	223K → 38K	1.86×

Form assembly is unchanged — identical kernel sizes and times to within noise.
gem.optimise.contraction already pulled these ListTensors out to be outermost so
that each entry could be factorised, so form compilation had already recovered the
structure by another route. The win is specific to dual evaluation.

Testing

Full FIAT suite. In Firedrake, matrix-free mult/multTranspose against the
assembled matrix, and an interpolate/restrict round trip, for NCF/NCE on hexahedra
and RTCF/RTCE on quadrilaterals, spectral and fdm variants.

Note on is_orthogonal

It now contracts the two maps with the gem interpreter, which neither belongs in an
element constructor nor scales. That is a deliberate temporary workaround and is
flagged with a TODO: the products of Deltas already say which components coincide,
so this should be read off symbolically, without the interpreter and without delta
elimination.

🤖 Generated with Claude Code

[pbrubeck]

Closed by mistake, not merged. Reordering the stack so that this PR sits below #268 put this branch's commit into pbrubeck/fix/dual-enriched, which was still this PR's base at the time; GitHub read that as merged and deleted the branch. A merged PR cannot be reopened.

Reopened as #271 — same branch, same commit (a75a2456), now based on pbrubeck/sum-factorise. Nothing here was merged into anything; please continue any review there.