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	#271	FInAT: select H(div)/H(curl) components with a Delta
3	#268	FInAT: dual evaluate on each sub-element's own points
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

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.

Supersedes #270, which GitHub marked merged by mistake:
reordering the stack put this branch's commit into pbrubeck/fix/dual-enriched,
which was still #270's base at the time, and a merged PR cannot be reopened.
Same branch, same commit, retargeted onto pbrubeck/sum-factorise.

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

[rckirby]

This seems reasonable. You haven't touched tests. Was this previously untested, or are there tests for correctness that still pass?

[pbrubeck]

This seems reasonable. You haven't touched tests. Was this previously untested, or are there tests for correctness that still pass?

I have the tests for these guys in #268