Facet QuadratureElement on TensorProductCell by pbrubeck · Pull Request #274 · firedrakeproject/fiat

pbrubeck · 2026-08-04T12:29:46Z

Interpolating into a Boundary Quadrature space on a quadrilateral raised

NotImplementedError: How to tabulate TensorProductElement on FacetPointSet?

A FacetPointSet maps the facet rule onto every facet at once. The facets of a
product cell lie in several directions, so those points have no single product
structure for a TensorProductElement to factor — factor_point_set cannot
split them.

The decomposition

The facets of a product lie in the direction of one factor at a time,

d(A x B) = dA x B  u  A x dB

so one direction's facets are a product even though all of them together are
not. QuadratureElement._facet_factorisation builds one
TensorProductElement per direction: the factor left out contributes its two
vertices — the whole of its own boundary — and the others contribute the
factors of the facet rule. Each summand's points are a TensorPointSet, which
is what dual evaluation needs of them, and the results stack along the basis
index.

The facet numbering already groups by direction, so the summands stack in
order:

UFCQuadrilateral   0 -> ((0,1),0)   1 -> ((0,1),1)   2 -> ((1,0),0)   3 -> ((1,0),1)
UFCHexahedron      0 -> ((0,1,1),0) ... 4 -> ((1,1,0),0)  5 -> ((1,1,0),1)

Why tabulation is left alone

The flat facet index is what lets the tabulation on each facet collapse into a
single Delta(facet_number, f) in _select_expression (the rule added in
#168). A direct sum splits that index per direction and loses it, so
tabulation keeps working on the FacetPointSet as it is.

Consequently the summands are not registered with as_enriched: they
number their degrees of freedom in each product's order, which differs from
this element's whenever the direction is not the first factor.

flat:      [f0q0 f0q1 | f1q0 f1q1 | f2q0 f2q1 | f3q0 f3q1]
sum dir 0:  index_shape (2, nq) -> f0q0 f0q1 f1q0 f1q1   agrees
sum dir 1:  index_shape (nq, 2) -> f2q0 f3q0 f2q1 f3q1   differs

dual_evaluation renumbers them by bringing the facet index outermost.
test_facet_factorisation_ordering pins that down for quadrilaterals and
hexahedra at several degrees.

Also here

A rule on a point is now a PointSingleton rather than a PointSet holding
one point, so the vertices a factor contributes carry no spurious unit index.
almost_equal moves up to AbstractPointSet to go with it.
TensorFiniteElement.dual_evaluation delegates to the base element on each
component instead of contracting the tensorised weights itself. The dual
basis is the base element's coupled to an identity on the value index, so it
is the same contraction — but only the base element knows the points its
functionals evaluate on, which is exactly what a facet quadrature element
needs to say. Without it a vector- or tensor-valued Boundary Quadrature
space on a quadrilateral fails the same way.

Status

Scalar, vector and tensor Boundary Quadrature interpolation now works on
quadrilaterals; values are exact at the points and agree from both sides of a
shared facet. Unchanged and still unsupported, as before this branch:
hexahedra (blocked by Firedrake's check_element) and genuine
TensorProductCells, where construct_subelement(int) raises because the
facets are not all the same shape.

🤖 Generated with Claude Code

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 · 2026-08-04T14:01:34Z

Paired firedrake PR: firedrakeproject/firedrake#5323 — firedrakeproject/firedrake#5323

It carries the tsfc zero propagation this needs to reach codegen, the regression tests for quadrilaterals, and a DROP BEFORE MERGE commit pinning CI to this branch.

pbrubeck · 2026-08-04T14:49:44Z

+        # Avoid circular import dependency
+        from finat.tensor_product import TensorProductElement
+
+        product = getattr(self.cell, "product", None)


Ugly getattr go away! Use isinstance(self.cell, FlattenedDimensions)

rckirby

Does this have top-level tests of correctness?

rckirby · 2026-08-05T20:39:31Z

        """GEM expression describing the points, with free indices
        ``self.indices`` and shape (point dimension,)."""

+    def almost_equal(self, other, tolerance=1e-12):


Do we want to read the default tolerance from some underlying data type to facilitate single precision in the future?

yes, that would be good. Currently we put a dtype on ref_el.vertices precisely becuase we only use tolerances for point inclusions/comparisons

rckirby · 2026-08-05T20:41:26Z

+
+
+@pytest.mark.parametrize("cellname", ["quadrilateral", "hexahedron"])
+@pytest.mark.parametrize("degree", [1, 2, 3, 5])


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>

Interpolating into a Boundary Quadrature space on a quadrilateral raised "How to tabulate TensorProductElement on FacetPointSet?". A FacetPointSet maps the facet rule onto every facet at once, and the facets of a product cell lie in several directions, so its points have no single product structure for a TensorProductElement to factor. The facets of a product lie in the direction of one factor at a time, d(A x B) = dA x B u A x dB, so one direction's facets are a product even though all of them together are not. Dual evaluate direction by direction: each is a product of quadrature elements whose points are a TensorPointSet, and the results stack along the basis index. Tabulation is left alone. It works on the FacetPointSet as it is, and the flat facet index is what lets the tabulation on each facet collapse into a single Delta(facet_number, f) in _select_expression; a direct sum would split that index per direction and lose it. So the summands here are not registered with as_enriched: they number their degrees of freedom in each product's order, which differs from this element's whenever the direction is not the first factor, and dual_evaluation renumbers them by bringing the facet index outermost. test_facet_factorisation_ordering pins that down. A rule on a point is now a PointSingleton rather than a PointSet of one point, so the vertices a factor contributes carry no spurious unit index; almost_equal moves up to AbstractPointSet to go with it. TensorFiniteElement.dual_evaluation now delegates to the base element on each component instead of contracting the tensorised weights itself. The dual basis is the base element's coupled to an identity on the value index, so this is the same contraction, but only the base element knows the points its functionals evaluate on, which is what a facet quadrature element needs to say. Without it a vector or tensor valued Boundary Quadrature space on a quadrilateral fails the same way. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

pbrubeck changed the title ~~Dual evaluate a facet quadrature element on a product cell~~ Facet quadrature element on TensorProductCell Aug 4, 2026

pbrubeck changed the title ~~Facet quadrature element on TensorProductCell~~ Facet QuadratureElement on TensorProductCell Aug 4, 2026

pbrubeck mentioned this pull request Aug 4, 2026

Interpolate into a Boundary Quadrature space on a quadrilateral firedrakeproject/firedrake#5323

Open

pbrubeck commented Aug 4, 2026

View reviewed changes

rckirby reviewed Aug 5, 2026

View reviewed changes

pbrubeck force-pushed the pbrubeck/tensor-facet-quadrature branch from 47c5baf to 975634c Compare August 18, 2026 10:03

pbrubeck force-pushed the pbrubeck/tensor-facet-quadrature branch 2 times, most recently from 7167c73 to a14c69c Compare August 24, 2026 21:15

pbrubeck force-pushed the pbrubeck/tensor-facet-quadrature branch from a14c69c to 121951c Compare August 26, 2026 11:11

pbrubeck and others added 9 commits August 26, 2026 12:49

style

4d52a34

calibrate tolerance

a099322

pbrubeck force-pushed the pbrubeck/tensor-facet-quadrature branch from 121951c to a099322 Compare August 26, 2026 11:57

Navigation Menu

Uh oh!

Facet QuadratureElement on TensorProductCell - #274

Facet QuadratureElement on TensorProductCell#274
pbrubeck wants to merge 9 commits into
pbrubeck/fix/dual-enrichedfrom
pbrubeck/tensor-facet-quadrature

pbrubeck commented Aug 4, 2026

Uh oh!

pbrubeck commented Aug 4, 2026

Uh oh!

pbrubeck Aug 4, 2026

Uh oh!

rckirby left a comment

Uh oh!

rckirby Aug 5, 2026

Uh oh!

pbrubeck Aug 7, 2026

Uh oh!

rckirby Aug 5, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants



		@pytest.mark.parametrize("cellname", ["quadrilateral", "hexahedron"])
		@pytest.mark.parametrize("degree", [1, 2, 3, 5])

Uh oh!

Conversation

pbrubeck commented Aug 4, 2026

The decomposition

Why tabulation is left alone

Also here

Status

Uh oh!

pbrubeck commented Aug 4, 2026

Uh oh!

pbrubeck Aug 4, 2026

Choose a reason for hiding this comment

Uh oh!

rckirby left a comment

Choose a reason for hiding this comment

Uh oh!

rckirby Aug 5, 2026

Choose a reason for hiding this comment

Uh oh!

pbrubeck Aug 7, 2026

Choose a reason for hiding this comment

Uh oh!

rckirby Aug 5, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants