Interpolate into a Boundary Quadrature space on a quadrilateral - #5323
Open
pbrubeck wants to merge 8 commits into
Open
Interpolate into a Boundary Quadrature space on a quadrilateral#5323pbrubeck wants to merge 8 commits into
pbrubeck wants to merge 8 commits into
Conversation
…ient compile_expression_dual_evaluation contracts the dual evaluation against the dual coefficient with an IndexSum over the basis indices when the dual argument is a Cofunction. When the target element has a block-structured dual basis (e.g. a RestrictedElement, whose dual evaluation is a gem.Concatenate), this hid the concatenation index from the subsequent unconcatenate() call, which can only split along indices that are still free in the assignment. Resolve the concatenation of the dual basis before the basis indices are summed over, following the same pattern already used for coefficient evaluation in tsfc.fem. This makes matrix-free adjoint interpolation (MatMultTranspose of a two-form Interpolate) work for restricted elements. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
EnrichedElement has no dual basis without it, so interpolating into the facet space of a hexahedron raises NotImplementedError, and gem cannot unconcatenate a dual basis whose blocks have already been contracted with a Cofunction. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
compile_expression_dual_evaluation carried a TODO to apply the GEM optimisations that assembly applies. Apply one of them: contract each assignment once unconcatenate has split the Concatenate nodes. This matters for H(div) and H(curl). Those elements are a sum of blocks, each mapping into one vector component, and with firedrakeproject/fiat#270 each block selects its component with a Delta. Before unconcatenation the blocks are branches of a Concatenate, so the Deltas never meet; afterwards each pair of blocks is its own assignment, the Deltas multiply, and delta elimination cancels every pair mapping to different components. Coarse to fine interpolation, 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 does not go through this function and is unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Apply to dual evaluation the zero handling assembly already does. Fold zero-valued Literals into symbolic Zero before anything else, so the algebraic identities built into gem (Product(_, Zero) -> Zero, and so on) can propagate them and the later steps can prune what they reach. Then drop any pair whose expression has collapsed to Zero once the blocks are contracted: a cross-component block that a Delta has just cancelled is a real zero contribution rather than one merely absent from the sum, and it must not reach codegen as a temporary. Emit a Zero temporary as a zero-filled read-only array when one does reach loopy, as gem.Zero carries no .array for the Constant branch to use.
Run the quadrature element tests on a quadrilateral mesh as well as a triangular one, and check that interpolation into a facet quadrature space reproduces the expression at the points and gives both cells sharing a facet the same values there. Scalar, vector and tensor valued spaces each take a different path through dual evaluation, so cover all three. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Point the existing FIAT install step at pbrubeck/tensor-facet-quadrature, which is pbrubeck/fix/dual-enriched plus the facet quadrature work. A Boundary Quadrature space on a quadrilateral cannot be interpolated into without it: the dual basis names the points of every facet at once, and a TensorProductElement cannot tabulate on those. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
pbrubeck
commented
Aug 4, 2026
pbrubeck
commented
Aug 4, 2026
pbrubeck
commented
Aug 4, 2026
Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>
pbrubeck
commented
Aug 4, 2026
pbrubeck
force-pushed
the
pbrubeck/hdivcurl-delta
branch
from
August 18, 2026 22:22
a992610 to
537d1b5
Compare
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
1 participant
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.
Paired with firedrakeproject/fiat#274 — firedrakeproject/fiat#274
FunctionSpace(quad_mesh, "Boundary Quadrature")could be assembled but notinterpolated into:
The element's dual basis names the points of every facet at once, and the
facets of a product cell lie in several directions, so those points have no
single product structure for a
TensorProductElementto factor. fiat#274 fixesthat by dual evaluating one facet direction at a time — each direction is a
product, since
d(A x B) = dA x B u A x dB.Here
tsfc: propagate zeros through dual evaluation— apply to dualevaluation the zero handling assembly already does. Fold zero-valued
Literals into symbolicZerofirst so gem's algebraic identities canpropagate them; drop pairs whose expression collapsed to
Zeroonce theblocks are contracted, since a
Delta-cancelled cross-component block is areal zero contribution rather than one merely absent from the sum; and emit a
Zerotemporary as a zero-filled read-only array if one does reach loopy,as
gem.Zerohas no.arrayfor theConstantbranch to use.Test Boundary Quadrature on quadrilaterals— the quadrature elementtests now run on a quadrilateral mesh as well as a triangular one, plus a new
interpolation test checking that the result reproduces the expression at the
points and that both cells sharing a facet see the same values there. Scalar,
vector and tensor valued spaces each take a different path through dual
evaluation, so all three are covered.
DROP BEFORE MERGE— points the existing FIAT install step atpbrubeck/tensor-facet-quadratureso CI builds against fiat#274.Verification
Against the base FIAT commit the three new
quadrilateralcases fail with theNotImplementedErrorabove; with fiat#274 the file is 25/25. Kernel cache wascleared for both runs — with a warm cache the failing cases pass spuriously.
Serial
test_interpolate,test_interpolate_zany,test_quadratureandtest_interpolate_vs_project: 106 passed, 0 failed. Parallel-marked tests werenot run locally.
Still unsupported, unchanged by this PR: hexahedra, which
check_elementrejects for this family, and extruded meshes, where the facets are not all the
same shape so the flat facet rule does not exist.
🤖 Generated with Claude Code