Interpolation into MixedSubElement - #261
Merged
Merged
Conversation
pbrubeck
force-pushed
the
pbrubeck/interp-mixed
branch
from
July 17, 2026 09:25
6d2573b to
25f3152
Compare
The transformers padded a ListTensor with zeros to place the value of the
wrapped element in the component it maps to. Contracting two of those
over the value index cannot cancel the components that do not coincide,
because the zeros are entries of a tensor rather than a property of the
expression, so every pair of blocks stays in the kernel.
Select the component with a Delta instead. 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 index.
Interpolating from NCF1/NCE1, with tsfc contracting the blocks once they
are separated, kernel temporaries and wall clock:
NCF 1->4 832K -> 68K 3.07x faster
NCE 1->4 307K -> 110K 2.23x faster
NCE 1->3 131K -> 61K 2.23x faster
Form assembly is unchanged; gem.optimise.contraction already pulled these
ListTensors out to be outermost so that each entry could be factorised.
is_orthogonal now contracts the two maps numerically, which neither
belongs in an element constructor nor scales; see the TODO there.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
find_group() collected indexed Concatenate nodes by intersecting against the free indices of the expressions, and asserted that the concatenation index was among them. But a Concatenate can only be split if the index is carried by an assignment *variable*: that is what split_variable() needs to slice the variable into blocks. Intersect against the variables' free indices instead. An index that has already been contracted away in the expression, or that belongs to another pair's variable, is now left alone rather than tripping the assertion or looping forever on a pair that cannot be split. The groups found are unchanged whenever the variables carry every concatenation index, which is the case for form compilation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Interpolating into the facet space of a hexahedron raised
"NotImplementedError: Dual basis not defined for element EnrichedElement",
and matrix-free adjoint interpolation out of a restricted space failed in
gem.unconcatenate. Both come from an element whose functionals do not all
evaluate on the same points having no way to be dual evaluated except
through dual_basis, which is a single weight tensor on a single point set.
Such an element is a direct sum. Give every element the sub_elements it is
the direct sum of -- itself, unless it is one -- and dual evaluate each on
its own points, stacking the results along the basis index the direct sum
occupies:
* an EnrichedElement is the sum of the elements it enriches;
* a TensorProductElement with a summed factor is the sum of the products
of that factor's sub-elements with the factors either side of it, since
the tensor product distributes over the direct sum;
* H(div)/H(curl), flattened and discontinuous wrappers carry the
decomposition of what they wrap.
The decomposition is one level deep and dual_evaluation recurses, so each
element concatenates along the index its own sum occupies. That matters for
a tensor product: its sub-elements own a slice of one factor's basis index
and the whole of every other, which is a sub-box of the basis multiindex
rather than a contiguous range of a flat one. Stacking on the factor's own
index is correct wherever the summed factor sits; flattening to leaves would
renumber the DoFs unless it were outermost.
Concatenating over the basis index is what unconcatenate is for: the index
stays free in the assignment, so it can be split downstream. A concatenation
over the points could not be -- points are always contracted by the sum in
dual_evaluation -- which is why the blocks are never brought onto common
points.
dual_point_set names the points of every sub-element, as a UnionPointSet, for
the callers that build a function space on them rather than contract against
them. It carries no structure across the sub-elements: contracting is the
business of each of them, and each keeps its own indices.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
EnrichedElement infers nodality from pairwise orthogonality, which only recognises blocks that map to different components. Restricting an element on a tensor product cell splits each direction into several blocks that map to the same component, so the inferred answer is False and the restriction of an NCE or NCF element was reported as non-nodal. Restriction selects disjoint subsets of the DoFs, so the restricted subelements are nodal whenever the original ones are, and a pullback applied blockwise preserves that. Carry the flag through instead of re-deriving it, as the tensor product restriction already does. The test exercises the decomposition several levels deep: the summed factor of these elements sits inside a FlattenedDimensions inside an H(curl) or H(div) wrapper, so reaching the elements that are not themselves a direct sum means recursing rather than reading one level. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
An element whose functionals do not all evaluate on the same points had no way to be dual evaluated: dual_basis is a single weight tensor on a single point set, and there is no such tensor for a direct sum. Interpolating into the facet space of a hexahedron raised "Dual basis not defined for element EnrichedElement", and interpolating into NCE of degree two or more failed in gem, as one factor of a summand is a direct sum in turn. Only a direct sum can evaluate summand by summand, so rewrite an element with a summed part as one, with the sum outermost. as_enriched() is that rewrite, and the sum commutes with everything that can hold it: * a pullback is linear, so it distributes over the sum; * a tensor product distributes over a sum in its first factor; * the wrappers that leave the dual basis alone are transparent. A tensor product and a pullback each delegate dual_evaluation to the rewrite, and EnrichedElement contracts each summand on its own points and concatenates along the basis index it occupies. That index stays free in the assignment, which is what unconcatenate needs to split it downstream; a concatenation over the points could not be split, as dual_evaluation contracts them away. Only the first factor of a product may be summed. Its summands own a contiguous range of the flat basis index and so stack in the order the product already numbers them, whereas a sum in a later factor would interleave with the factors before it and stacking would renumber the degrees of freedom. Every element of these families is built with the summed factor first. EnrichedElement now has a dual_basis, block diagonal on the union of its summands' points, for the callers that want a single weight tensor rather than to contract against one. It blocks over the same summands dual_evaluation does, so the two cannot disagree about which points each summand owns. A UnionPointSet is a PointSet that records where each summand's points begin, so a function space can be built on it. An element that needs structure to tabulate -- a product, which cannot factor a union -- tabulates on each point set of the union in turn rather than on the union unstructured, keeping the structure each of them has to be sum factorised. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A QuadratureElement whose rule is a union of point sets has no single set of basis functions evaluating on all of them, so each summand evaluates on its own points and the results stack along the basis index. Tabulation matches `ps` against the rule's blocks in the same way, and a vector or tensor wrapper of a direct sum distributes over the sum first. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
FiniteElementBase.dual_evaluation brings the sum outermost and delegates to a _dual_evaluation hook, so as_enriched is called from one place rather than from each element that can hold a sum. Drop the overrides in TensorProductElement and WrapperElementBase, and reduce the ones in TensorFiniteElement and EnrichedElement to the hook. Fix QuadratureElement's union branch, which unpacked a pair from the triple dual_evaluation returns, and returned a pair in turn. Firedrake builds a quadrature space on the points of the target's dual basis, so this is reached whenever that target is a direct sum; test it. Collapse restrict_hcurl and restrict_hdiv onto WrapperElementBase, and trim the docstrings to :arg: style. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>
A `QuadratureElement` whose rule sits on a `UnionPointSet` split the union
itself, in its own `_dual_evaluation`. That only covered the element as it
stands: a `TensorFiniteElement` wrapping it takes the weights from
`dual_basis` and calls the callable once, on the whole union. The callable
then tabulates into a `Concatenate` over the points, which is contracted
away before `unconcatenate` can split along it, and the `Concatenate`
reaches loopy:
AssertionError: cannot generate expression from <class 'gem.gem.Concatenate'>
Firedrake interpolates into a vector-valued quadrature space on the points
of the target's dual basis, so every H(div)/H(curl) grid transfer took that
path.
Register the union with `as_enriched` instead. `dual_evaluation` brings a
direct sum outermost before it evaluates, so the summands then each get
their own point set wherever the element appears, the tensor wrapper
included, and the hand-rolled split goes away.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
pbrubeck
force-pushed
the
pbrubeck/interp-mixed
branch
from
August 27, 2026 15:40
ddc4c48 to
564a39d
Compare
pbrubeck
force-pushed
the
pbrubeck/interp-mixed
branch
from
August 27, 2026 15:44
564a39d to
bafa803
Compare
A tensor element couples a basis function to a value component through an identity, so it needs no rewrite of its own to reach a direct sum: dual evaluating the base element one component at a time gives the components whatever structure the base element has, a direct sum among it, and keeps the basis indices in the element's own index shape rather than the flat one a sum of tensor elements gives. The identity is contracted rather than carried into the basis indices, as a free index into a ListTensor defeats argument factorisation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
pbrubeck
force-pushed
the
pbrubeck/interp-mixed
branch
from
August 27, 2026 15:49
bafa803 to
b01a523
Compare
Fix the getattr/setattr lazy init in QuadratureElement._weights, dedupe the as_enriched summand-distribution shape into distribute_over_sum, rename find_group's free_indices to splittable_indices, replace the repeated isinstance(ps, UnionPointSet) checks with a point_sets property, restore DiscontinuousElement's dual evaluation as the _dual_evaluation hook instead of a public override, and dedupe repeated test setup. Guard EnrichedElement.dual_basis to raise NotImplementedError for a summand with its own dual_transformation, since the property has no coordinate_mapping to apply it with and would otherwise silently combine untransformed weights. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BjZ6kYK1tmLfm6wUxXwVtG
A sum whose summands are not mutually nodal has no dual basis to give. Its functionals are still the block diagonal weights that dual_basis builds, but every caller contracts those into coefficients, and for a non-nodal sum that contraction is the inverse of the nodal matrix, not the identity. MINI is the case in hand: the linear basis functions are not zero at the bubble's point. dual_evaluation already refused this. dual_basis did not, so the sum advertised weights that no caller could use, and has_pointwise_dual_basis answered False only because the block diagonal carries Zero blocks rather than because the sum is non-nodal. Refuse it in dual_basis as well. dual_evaluation keeps its own guard: it contracts each summand through that summand's dual_basis and never asks self for one. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013NoTXWyj2fVdJTHnMDFB4k
is_orthogonal only special-cases HCurlElement/HDivElement pairs, so it conservatively returned False for MixedSubElement summands and blocked EnrichedElement._dual_evaluation with NotImplementedError. A MixedSubElement embeds each summand into its own disjoint slot of the flattened value shape, so they are always orthogonal by construction. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01RFoiWhAUbP53byrgoqxx7H
Phase 3 of _collect_monomials summed the indices that no atomic carries without cancelling the Deltas over them first, and left the remaining factors as the product trees a nested expansion had already built, where a Delta is invisible. Such a contraction was carried out rather than cancelled, keeping an index that the Delta should have renamed away. Flatten the non-atomic factors and eliminate Deltas over every index the monomial contracts, as gem.optimise.contraction already does for coefficient evaluation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015Vnd4VHtANbfmiMxSPSXSG
Phase 3 of _collect_monomials flattens the non-atomic factors before cancelling Deltas, and did so in two unsound ways. traverse_product started a fresh renamer, so a contraction it hoisted out of the rest kept an index the monomial was already summing over, and the duplicate tripped MonomialSum.add's assertion. Hand it the renamer that has already seen the monomial's own indices. The product being flattened also folds to Zero as soon as one of its factors is Zero, taking the free indices of the others with it, and sum_factorise was then left to contract an index that no factor carried. Such a monomial is zero, so recognise it and skip the contraction rather than planning one. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01F1K5upgeZTKWj9m4tHt9Vn
Phase 3 of _collect_monomials flattens the non-atomic factors before cancelling Deltas, and did so in two unsound ways. traverse_product descends into contractions, so flattening pulled apart the coefficient evaluations that gem.optimise.contraction had already planned and re-planned them as a single problem. An interior facet integral of two DG3 coefficients on a hexahedron then handed sum_factorise one connected contraction over eight indices -- the two facet quadrature indices and each cell's three basis indices -- which it refuses. Stop the flattening at the contractions the product already contains; the Deltas that have to be cancelled are factors of the product tree, not of those. The product being flattened also folds to Zero as soon as one of its factors is Zero, taking the free indices of the others with it, and sum_factorise was then left to contract an index that no factor carried. Such a monomial is zero, so recognise it and skip the contraction rather than planning one. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01F1K5upgeZTKWj9m4tHt9Vn
collect_monomials cancels Deltas over the indices contracted by a monomial, which needs the non-atomic factors flattened first so that the Deltas expansion buried inside a product are visible. Rebuilding rest from that flattened list unconditionally discards the grouping the recursion had already chosen, and sum_factorise re-associates it worse: hexahedral Q3 Helmholtz cost 2.3% more flops, and every laplacian, Helmholtz and H(div) kernel measured paid between 0.2% and 2.3%. Cancel in _cancel_deltas, which hands back rebound factors only when a Delta does cancel, and otherwise the factors it was given. Form assembly flop counts then match what they were before the cancellation was introduced, and interpolation against a Cofunction dual argument keeps its win: 1128 -> 360 flops on a degree 4 quadrilateral, 140039 -> 25000 on a degree 8 facet-restricted hexahedron, with no Delta left for the replace_delta fallback to lower into an identity matrix. Take the sum indices that reach the rest from the free indices of its factors, rather than as those the atomics do not carry. A Zero absorbs the free indices of the factors it multiplies, so the two differ exactly when the monomial is zero, and sum_factorise then contracts nothing and folds the Zero on its own. Drop the index_replacer argument to traverse_product, which stop_at makes unreachable, and pin the cancellation with a test: main carries the Delta into the contraction, and flattening alone leaves the nested one behind. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0158RTiFnSyhgKbsbz11JMak
Both sides rewrote how phase 3 of _collect_monomials flattens the non-atomic factors before cancelling Deltas. Keep this side's flattening, which stops at the contractions the product already contains. The other side descends into them and threads a renamer through traverse_product so that the indices it hoists stay distinct, but descending re-plans the contractions gem.optimise.contraction had already planned, and sum_factorise re-associates them worse: hexahedral Q3 Helmholtz costs 2.3% more flops, and every laplacian, Helmholtz and H(div) kernel measured pays between 0.2% and 2.3%. Nothing is hoisted when the flattening stops, so the indices cannot collide and the renamer has nothing to keep apart. Take the pruning of the cancelled factors, which delta_elimination leaves wrapped in a ComponentTensor that nothing downstream is obliged to remove. Two survived per dual evaluation kernel here; none now, at no cost in flops. The Zero case both sides special-cased is already covered: this side takes the sum indices that reach the rest from the free indices of its factors, so a Zero that has absorbed them contracts nothing and folds on its own. Keep the traverse_product signature from the other side, whose renamer argument subsumes the rename_map it replaces. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0158RTiFnSyhgKbsbz11JMak
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01G4zn3Z5RoNPnaVpDzqP2dg
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01G4zn3Z5RoNPnaVpDzqP2dg
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01G4zn3Z5RoNPnaVpDzqP2dg
pbrubeck
added a commit
to firedrakeproject/firedrake
that referenced
this pull request
Sep 9, 2026
This PR needs firedrakeproject/fiat#261 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. The step runs after Firedrake is installed, so that resolving its dependencies cannot pull fenics-ufl back from FEniCS/ufl@main over the branch, and calls firedrake-clean so no kernel cached against the released UFL survives. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
pbrubeck
added a commit
to firedrakeproject/firedrake
that referenced
this pull request
Sep 9, 2026
This PR needs firedrakeproject/fiat#261 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. The step runs after Firedrake is installed, so that resolving its dependencies cannot pull fenics-ufl back from FEniCS/ufl@main over the branch, and calls firedrake-clean so no kernel cached against the released UFL survives. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
pbrubeck
added a commit
to firedrakeproject/firedrake
that referenced
this pull request
Sep 9, 2026
This PR needs firedrakeproject/fiat#261 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 UFL. A step of its own after the install therefore never gets to run. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
An error occurred while trying to automatically change base from
pbrubeck/fix/dual-enriched
to
pbrubeck/hdivcurl-delta
September 10, 2026 11:29
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Uh oh!
There was an error while loading. Please reload this page.