Sum factorise independent contractions separately

[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.

Interpolating a product of three coefficients into a hexahedron CG1 space:

mesh = ExtrudedMesh(UnitSquareMesh(1, 1, quadrilateral=True), 1)
V = FunctionSpace(mesh, "CG", 1)
Function(V).interpolate(Function(V)*Function(V)*Function(V))

raised

File "gem/optimise.py", line 397, in sum_factorise
    raise NotImplementedError("Too many indices for sum factorisation!")

sum_factorise searches every ordering of the contraction indices, which is
factorial in their number, so it gave up past six of them. On a tensor product cell
each coefficient evaluation contributes one index per direction, so three
coefficients on a hexahedron already need nine.

But those indices do not interact: no factor carries the indices of more than one
coefficient, so the orderings that interleave them are never worth searching. This
splits the contraction into the connected components of the graph joining indices
that share a factor, and searches each component separately.

That fixes the failure above, and is also cheaper wherever it already worked, since
the search is over the orderings of each component rather than of all the indices at
once. A single connected contraction is still bounded at six.

Testing

test/gem/test_sum_factorise.py checks the factorised expression against the value
obtained by contracting the tables directly, for products of up to five independent
contractions of up to five indices each — well beyond the old limit — and checks that
one connected contraction is still rejected. Both new cases were confirmed to fail on
the unpatched code with the original error.

Verified end to end in Firedrake that f**n interpolates correctly for n up to 5 on
quadrilateral, triangle, hexahedron and prism meshes, to machine precision.

Full FIAT suite: 2431 passed, 26 skipped, 31 xfailed.

🤖 Generated with Claude Code

[pbrubeck]

@connorjward this is a bugfix, should I target release?

[pbrubeck]

I think that we only hit this in interpolate and not in assemble because assemble does some extra COFFEE processing. But arguably the right thing is to have the decomposition into smaller subproblems in one place.

[connorjward]

@connorjward this is a bugfix, should I target release?

I don't especially mind. Given that we are releasing 2026.10.0 soonish, and I have no intention of making another patch release unless required, the issue is somewhat moot.

[rckirby]

I think this is fine, if kind of complicated. Is there a test that this is correct, or at least returns gem that evaluates to something equivalent to the previous version?