Sum-factorise interpolation into composite EnrichedElements by pbrubeck · Pull Request #294 · firedrakeproject/fiat

pbrubeck · 2026-09-10T15:05:18Z

Interpolation into facet-restricted finite element spaces on tensor-product cells was not being sum-factorised correctly. Three separate problems contributed.

1. Nested direct sums were not fully expanded

An enriched element can contain another direct sum. This occurs, for example, in MixedElement(EnrichedElement(...)). Different code paths used different decompositions of the element. Basis tabulation stopped at the immediate elements, while dual evaluation could expand nested sums. The basis and dual basis then had different blocks, so TSFC could not pair the corresponding pieces.

The element now exposes one fully expanded list of direct-sum leaves. Basis tabulation, point evaluation, dual evaluation, and dual-basis construction all use that list and the same ordering. Wrappers such as FlattenedDimensions are distributed over the sum so each leaf remains on the correct cell.

2. All-zero tables were not folded

_constant_fold_zero_literal used numpy.array_equal(array, 0). For a non-scalar array, that compares the array shape with the scalar shape and returns false. Dense tables containing only zeros therefore survived optimisation.

These tables are common for direct sums. Each component is zero at the points belonging to the other components. The pass now checks the values directly and preserves the literal dtype. Later simplification can then remove the zero Indexed and Product expressions.

3. Direct-sum contractions were not split

In simple terms, Concatenate(A, B) is one long vector made by putting A and B next to each other. If two such vectors are multiplied entry by entry and summed, matching blocks can be summed separately:

sum_j [A; B]_j [C; D]_j
  = sum_ja A_ja C_ja + sum_jb B_jb D_jb

There are no cross-terms because each position in the first vector meets only the same position in the second vector. split_contraction applies this identity to the direct-sum blocks produced by FInAT. It lets TSFC contract each block over its own points and recursively handles more than one summed index.

EnrichedElement._dual_evaluation now leaves the summand point indices free. TSFC can therefore apply split_contraction and choose the correct contraction for each block.

Validation

The FInAT dual-basis tests cover nested direct sums and preserve the summand point indices.
A regression test covers dual evaluation of an extruded-cell NCE x Q mixed element, where NCE is a composite EnrichedElement.
GEM tests cover folding zero tables, retaining nonzero tables, and removing products that contain zero tables.
Firedrake TSFC tests cover sum-factorised interpolation into facet-restricted spaces on tensor-product cells.

AI tools used in preparing this change: OpenAI Codex and Claude Code (Opus 5).

This PR needs firedrakeproject/fiat#294 and the UFL work in FEniCS/ufl#511 and FEniCS/ufl#512, none of which has merged, so CI installs both branches over the pins in pyproject.toml. Drop this commit once they land; nothing else on the branch touches .github. Both installs sit inside the Install Firedrake step, ahead of the Firedrake install rather than after it: that step ends with firedrake-clean, which imports Firedrake, and Firedrake cannot be imported against the released dependencies. A step of its own after the install therefore never gets to run.

numpy.array_equal compares shapes before it compares values, so it returned False for every table that the pass gave it and folded nothing but a scalar zero. A table of zeros stayed a dense Literal. Dual evaluation between facet-restricted elements on a tensor product cell tabulates the interior basis functions of each direct-sum component at the boundary points of every other component. Those tables are exactly zero, so 30 of the 49 component pairs that a hexahedron produces contribute nothing. Folding the tables lets Indexed and Product drop those pairs, which returns the interpolation to the O(degree^(dim + 1)) cost that sum factorisation gives. The kernel of tests/tsfc/test_dual_evaluation.py::test_dual_argument_is_sum_factorised loses 72% of its flops at degree 16. Fold on the value of the table, and keep the dtype so that a table of integers does not become a floating point zero. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013tWTW5ErfhCgV63nXai8HN

EnrichedElement presented two decompositions of the same direct sum. basis_evaluation concatenated blocks of shape elem.index_shape over self.elements, giving (4,4,3), (24,4) for NCE3, while _dual_evaluation concatenated whatever the as_enriched rewriting returned, giving (4,4,3), (96,). A basis and its dual basis must be blocked alike, so nothing downstream could pair them up and contract them. Promote EnrichedElement's private _summands to a `summands` property on every element, and block basis_evaluation, point_evaluation, dual_basis and _dual_evaluation along it alike. as_enriched on a FlattenedDimensions dropped the wrapper and returned summands on the tensor product cell, which cannot tabulate against the entities of the quadrilateral or hexahedron they came from. Distribute the flattening over the sum instead, as the other wrappers already do. Add split_contraction, which carries the identity that a sum over a whole direct sum is the sum of the sums over its blocks. Unlike unconcatenate it needs no assignment variable to carry the concatenation index, because the sum itself is what the Concatenate splits against. split_group holds the part that the two now share. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01KWajdMc5VFPbPuB1cupu9F

rckirby

A lot of the verbiage in commentary is vague/awkward. Please clarify.
Also, I see some new tests. Can we confirm that the feature being added or that was previously broken is now being tested?

rckirby · 2026-09-11T18:07:58Z

-        downstream; a concatenation over the contracted points could not be.
+        The summands do not share their points, so their evaluations stack
+        along the basis index while retaining their own point indices.
+        Concatenating over a free basis index is what


"Concatenating" is what "unconcatenate" does is confusing to me.

rckirby · 2026-09-15T14:15:49Z

-        summand may be a direct sum in turn.  These are the elements that
-        evaluate their dual basis on their own points, and whose points make
-        up the union that :attr:`dual_basis` works against.
+        An element is brought out as a direct sum one level at a time.  A summand


One level at a time? Clarify please?

The fully-flattened decomposition of this element as a sum of non-EnrichedElements

rckirby · 2026-09-15T14:19:31Z


+    @cached_property
+    def summands(self):
+        """The direct summands whose bases stack into this element's basis.


Clearer: Return (E1, E2, ..., En) where E = E1 \oplus E2 \oplus ... \oplus En?
Is this the correct idea?

rckirby · 2026-09-15T14:25:18Z

+def _unconcatenate(cache, pairs):
+    # Tail-call recursive core of unconcatenate.
+    # Assumes that input has already been sanitised.
+    # Only an index carried by an assignment variable can be split against it.


Can this be stated more clearly?
What happens if this fails? Is there a check?

rckirby · 2026-09-15T14:25:45Z

+    index, multiindices, mappings = split_group(cache, concat_group)
+
+    def cut(node):
+        """No need to rebuild expression of independent of the


Grammar? "expression of independent of the ..."

this grammar mistake was carried over

rckirby · 2026-09-15T14:26:59Z

+def split_contraction(expression, indices, cache=None):
+    """Splits a contraction along the :py:class:`Concatenate` nodes it sums over.
+
+    No assignment variable need carry the concatenation index here.  The sum


This is akward wording. Clarify, the math below is helpful!

Preserve direct-sum point indices in dual evaluation

808ca7e

pbrubeck added the LLM used An LLM was used in the production of this PR label Sep 10, 2026

pbrubeck mentioned this pull request Sep 10, 2026

TSFC: support same-mesh ufl.Interpolate in Form kernels firedrakeproject/firedrake#5258

Open

pbrubeck and others added 2 commits September 11, 2026 23:40

rckirby requested changes Sep 15, 2026

View reviewed changes

pbrubeck added 3 commits September 15, 2026 16:16

some review comments

e76f22b

Clarify direct-sum review comments

a36f371

Trim direct-sum docstrings

e0fa4fa

pbrubeck requested a review from rckirby September 15, 2026 16:01

Add direct-sum contraction regression tests

ff288ff

pbrubeck changed the title ~~Preserve direct-sum point indices in dual evaluation~~ Sum-factorise interpolation into composite EnrichedElements Sep 15, 2026

Cover NCE Q mixed dual evaluation

143ded3

pbrubeck commented Sep 15, 2026

View reviewed changes

Comment thread test/finat/test_dual_basis.py Outdated

pbrubeck added 2 commits September 15, 2026 17:26

Apply suggestion from @pbrubeck

906e2ca

Assert direct-sum point indices stay free

ff6e132

Navigation Menu

Sum-factorise interpolation into composite EnrichedElements - #294

Sum-factorise interpolation into composite EnrichedElements#294
pbrubeck wants to merge 10 commits into
mainfrom
pbrubeck/fix/dual-enriched-again

pbrubeck commented Sep 10, 2026 •

edited

Loading

Uh oh!

rckirby left a comment

Uh oh!

rckirby Sep 11, 2026

Uh oh!

Uh oh!

rckirby Sep 15, 2026

Uh oh!

pbrubeck Sep 15, 2026

Uh oh!

Uh oh!

rckirby Sep 15, 2026

Uh oh!

pbrubeck Sep 15, 2026

Uh oh!

rckirby Sep 15, 2026

Uh oh!

rckirby Sep 15, 2026

Uh oh!

pbrubeck Sep 15, 2026

Uh oh!

pbrubeck Sep 15, 2026

Uh oh!

rckirby Sep 15, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

pbrubeck commented Sep 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

1. Nested direct sums were not fully expanded

2. All-zero tables were not folded

3. Direct-sum contractions were not split

Validation

Uh oh!

rckirby left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

pbrubeck commented Sep 10, 2026 •

edited

Loading