Automate zany basis transformations by pbrubeck · Pull Request #259 · firedrakeproject/fiat

pbrubeck · 2026-07-13T09:36:12Z

Add finat/zany.py, to automatically derive the basis transformation V = E V^c D of Kirby (2017) directly from a FIAT element's dual basis:

Automated

AI tool: Claude Fable 5

Add finat/zany.py, a prototype that derives the basis transformation V = E V^c D of Kirby (2017) directly from a FIAT element's dual basis: - Sparsity comes from entity_dofs() and functional types: push-forward invariant nodes (point evaluations, integral averages) get identity rows. - The chain-rule factor V^c for facet normal-derivative moments reduces to Gram algebra, a = detJ sqrt(det Ghat / det G), b = G^{-1} T^T J n, in the frame of the physical normal and mapped reference tangents, valid in any dimension with no orientation logic. - The interpolation factor D is computed numerically by dual evaluation of tangential completion functionals (built from the normal node's own quadrature rule) against the nodal basis, replacing hand-derived univariate exactness rules. - The conditioning h-scaling generalizes via each node's max_deriv_order. The automated Morley transformation matches the hand-coded element to machine precision in 2D and 3D with a single dimension-independent code path, and passes check_zany_mapping. Move check_zany_mapping/make_unisolvent_points into the finat conftest and provide the checker as a fixture, since --import-mode=importlib forbids imports between test modules. Also add AGENTS.md notes on the transformation theory and the implementation conventions. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Replace the per-functional-kind row builders with a generic framework: - finat/functional.py: a degree of freedom is represented symbolically as l(f) = sum_q w_q <D, grad^m f(x_q)>, constructed from any FIAT functional's pt_dict/deriv_dict with the direction recovered numerically, so no dispatch over FIAT functional types is needed. Operations: covariant pullback (direction contracts with J), numeric dual evaluation against a nodal basis, and direction replacement. - finat/zany.py: FacetFrame holds the normal/tangential frame and a generic adjugate-based symbolic solve; one assembly loop treats every node identically. Push-forward invariance of value and tangential nodes is derived, not type-checked. Completion remainders recurse through already-assembled rows of V, ordered by entity dimension. - finat/morley.py: Morley is reimplemented on the framework; its basis transformation is one call to zany_basis_transformation. morley_transform moves verbatim to finat/walkington.py, its only remaining user. The remainder coefficients r_k = x_k - c beta_k are invariant under the sign ambiguity of the numerically recovered direction; the naive a*x_k formula flips sign exactly on the edges where the SVD picks the opposite orientation. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Derivative nodes away from facets have no geometric frame: FIAT keeps Cartesian directions on the physical cell, so the group of derivative nodes on an entity acts as its own completion (this is exactly the affine-interpolation equivalent case). The pulled-back direction is expanded in the group's own numeric direction basis, with weight-ratio factors keeping the expansion invariant under the scale and sign ambiguity of each node's recovered direction factorization. finat.Hermite now derives its transformation automatically; the matrix matches the deleted hand-coded one exactly in 2D and 3D, scaled and unscaled. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

pbrubeck · 2026-07-13T10:12:27Z

+        row[i] = c
+        for k, that in enumerate(frame.tangents):
+            r = x[k + 1] - c * beta[k]
+            coefficients = ell.with_direction(that).evaluate(fiat_element)


We should avoid evaluating the entire set of basis functions. We should only evaluate the basis subset suported on this facet.

Might not be obvious -- it may vary based on what we're evaluating on a facet (e.g. values vs derivatives)

That is true. In that case we should collect all points we plan to tabulate on and only call tabulate once as we do in toriesz()

Yes, and plus, fiat's tabulation is still going to shell out to the ON basis and Vandermonde conversion.

I think the savings of finding sparsity here versus the extra code complexity, but I could be wrong.

pbrubeck · 2026-07-13T10:48:04Z

+"""Symbolic representation of degrees of freedom.
+
+A :class:`Functional` represents a degree of freedom in the form
+
+.. math:: \\ell(f) = \\sum_q w_q \\langle D, \\nabla^m f(x_q) \\rangle,


This was a human suggestion midway through the planning/design stage. AI had already been considering it. I had to explicitly drive the agent to take this path.

pbrubeck · 2026-07-13T10:50:13Z

+class Functional:
+    """Symbolic degree of freedom with a single derivative direction.


The symbolic finat.Functional could potentially remove the inverse(M.T) in PhysicallyMapped.dual_transformation

Not sure what this means? the dual transformation acts on collections of such functionals.

Generalize Functional to arbitrary derivative order: directions live in derivative multi-index space, and pullback distributes them over a symmetric tensor, contracts each slot with the Jacobian, and collapses back. Point-jet groups are split per order so that Argyris/Bell vertex gradients and Hessians each solve in their own direction basis. Facet completions now recurse through vertex-jet rows and same-edge trace moments with no new code. zany_basis_transformation gains two options: avg=False divides facet-moment columns by the physical facet measure (the legacy convention where physical edge moments are plain integrals), and ndof drops trailing columns so constrained elements (Bell) can discard the constraint dofs of their extended element. The automatic matrices match the deleted hand-coded ones to machine precision for Argyris degrees 5-7 (integral variant with both avg conventions, and the point variant) and Bell. One convention change: the generic h^-m conditioning scaling now also applies to integral Argyris edge moments, which the hand-written code left unscaled; this is invisible when cell_size == 1. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

…rcier Extend the symbolic Functional with a value rank, parsing component weight profiles from pt_dict. Under the contravariant Piola map the scalar roles are mirrored: scaled facet normals are cofactor images K = adj(J)^T of the reference ones (exactly the physical compute_scaled_normal), so pure normal moments are invariant, while scaled tangents map by J. _piola_facet_rows matches per-point frame-coordinate profiles within each facet group (covering 3D MTW's point-varying RT-mapped tangential directions) and eliminates the residual normal profile through the Vandermonde recursion. FIAT builds 3D tangential components on the reciprocal basis cross(n, t_k), which transforms in-plane contravariantly; the tangent Gram change S is absorbed into the coordinate mixing (S = 1 in 2D). Interior moments are Piola-invariant by construction. Double contravariant tensor elements reuse the same code with one contraction per value slot. MardalTaiWinther (2D/3D, orders 1-2) and JohnsonMercier (2D/3D) now derive their transformations automatically, matching the deleted hand code to machine precision. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Vertex and edge point values of contravariant Piola elements are value point-groups whose components pull back through the cofactor matrix K = adj(J)^T, mirroring the Cartesian point-jet case. The trailing tangential facet constraints of the extended element are dropped via ndof, reusing PiolaBubbleElement's reduced entity_dofs bookkeeping. The hand-derived vertex-facet coupling correction of PiolaBubbleElement emerges automatically from the Vandermonde residual elimination. The automatic matrices match the hand-coded ones to machine precision for orders 0-2 in 2D and 3D. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

…xins Replace the free function zany_basis_transformation and its internal "if piola" branch with a template-method design: - PhysicallyMappedElement.basis_transformation (physically_mapped.py) now implements the entity-by-entity assembly loop directly, calling four hooks -- _check_mapping, _invariant_dofs, _facet_dof_rows, _point_dof_rows -- that carry all mapping-specific knowledge. The loop itself has no branch on the kind of pullback. - finat/zany.py supplies the two mixins implementing those hooks: ScalarPhysicallyMappedElement (affine pullback) and PiolaPhysicallyMappedElement ((double) contravariant Piola), plus the pure math functions they call (FacetFrame, _scalar_facet_rows, _scalar_point_rows, _piola_facet_rows, _piola_point_rows), which take plain arrays/GEM expressions with no self, keeping the mathematics readable independent of the class plumbing. - Concrete elements (Morley, Hermite, Argyris, Bell, MardalTaiWinther, JohnsonMercier, GuzmanNeilanFirstKindH1) drop their basis_transformation override entirely: mixing in the right base class is enough. The ndof truncation parameter is gone; the loop always slices by self.space_dimension(), which constrained elements already override. - finat.Functional is renamed to finat.PhysicallyMappedFunctional. Full finat suite, flake8, and pydocstyle all pass unchanged; the automatic matrices are byte-for-byte the same computation as before, just reorganized. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

PhysicallyMappedElement (physically_mapped.py) is restored to its original shape: a generic, abstract mixin with no knowledge of the zany theory, exactly as used by hand-coded elements (AW, HCT, PowellSabin, Walkington, ...). It no longer imports finat.functional. The template-method loop and its four hooks move into a new ZanyPhysicallyMappedElement(PhysicallyMappedElement) in finat/zany.py, which ScalarPhysicallyMappedElement and PiolaPhysicallyMappedElement now subclass instead of PhysicallyMappedElement directly. This keeps all zany-specific machinery -- the loop, the hooks, the pure math functions, and the symbolic PhysicallyMappedFunctional dependency -- contained in finat/zany.py and finat/functional.py, with physically_mapped.py staying general-purpose infrastructure. No behavioral change: full finat suite, flake8, and pydocstyle pass unchanged. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Replace the pre-project rough draft with concrete, grounded strategies: - Mathematical structures: the five-number dof shape independent of FIAT class, pullback as a uniform tensor-slot contraction (J for derivatives, K = adj(J)^T for Piola values), frame decomposition as the one shared computational primitive, extended elements as restriction. - Confusing mathematical ideas clarified: the papers' vs FInAT's Jacobian direction, the reciprocal-basis contragredient transformation that only shows up in 3D Piola elements (not in any of the papers), and the sign/scale invariance required of formulas built from numerically-recovered (SVD) invariants. - Confusing GEM patterns clarified: operator overloading vs manual node construction, adjugate/determinant as the symbolic matrix inverse, numpy.full with Zero() for sparse GEM arrays, the identity-then-mutate assembly idiom, and the row/column/transpose convention. - Design strategies that generalize: duality-first design, recovering mathematical type from data instead of class hierarchies, template method for mapping-specific case-splits, validating against independently computed ground truth, and incremental generalization from one worked example at a time. - A dedicated subsection on working with the human collaborator: the literature is necessary but not sufficient (the reciprocal-basis fix came from the user, not the papers), expecting and re-verifying after mid-project architectural corrections, and keeping the *why* on record for future sessions. Also fixes a stale pointer: morley_transform now lives in finat/walkington.py, not finat/morley.py (moved in an earlier commit when Morley itself was automated). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

rckirby

This is taking some work to get one's head around, but seems like an excellent initial stab. I'm still trying to digest what exactly going on, but passing tests are a great sign.

rckirby · 2026-07-13T15:48:49Z

+        return type(self)(self.points, self.weights,
+                          order=self.order, direction=direction)
+
+    def pullback(self, J: Node) -> "Functional":


Math for what this does? If we are acting on a physical-space function by changing its coordinates, then don't we need the whole coordinate mapping and not just its Jacobian?

Also, a comment that this does what it does whether or not the element is being mapped by coordinate change/Piola/etc. could be in order?

Finally, can this work with spatially varying or only constant coefficient Jacobians? (worth documenting...)

rckirby · 2026-07-13T15:58:59Z

+        -------
+        Functional
+            The functional with direction :math:`J \\otimes \\dots
+            \\otimes J : D`, acting on physical derivatives at the


Is this transforming the functional or derivatives? Please smooth out the documentation!

D defines the directional derivate. I think to answer your question you need to fix a frame of reference.

rckirby · 2026-07-13T15:59:43Z

+            direction[k] = direction[k] + T[index]
+        return self.with_direction(direction)
+
+    def evaluate(self, fiat_element: FiniteElement) -> numpy.ndarray:


apply_to_basis might be a better name? evaluate sounds like it should eat a function and spit out a number.

rckirby · 2026-07-13T16:00:28Z

-            cur += d
-
-        return ListTensor(M)
+    def basis_transformation(self, coordinate_mapping: PhysicalGeometry) -> ListTensor:


I approve of this massive deletion!

rckirby · 2026-07-13T18:41:46Z

+  already-assembled rows of $V$ (entities processed in increasing dimension), which
+  will later let completions couple to vertex jets (Argyris/HCT) for free.
+
+Derivative nodes *away from facets* (`_scalar_point_rows`, covering Hermite vertex


"away from facets"? Meaning on a vertex (which is on multiple facets...)?

rckirby · 2026-07-13T18:42:23Z

+
+Key facts the framework rests on:
+
+* **FIAT normals are "UFC consistent":** computed from the tangents by the same formula


FIAT normals are only UFC consistent if you use one of the UFC reference cells. Other cells are possible. Does that mean the theory breaks in that case?

rckirby · 2026-07-13T18:43:14Z

+  (assumes $\det J > 0$). For the record, the fully simplified closed forms the solve
+  reproduces are $a = \det J\sqrt{\det\hat G/\det G}$ and $b = G^{-1}T^TJ\hat n$ with
+  Gram matrices of the (mapped) tangents.
+* **Integral averages are push-forward invariant.** FIAT's Morley dofs are averages


Without type information, how are we determining that a given linear combination of point values is a scaled or unscaled moment?

rckirby · 2026-07-13T18:48:10Z

+   `entity_dofs()` and functional types (`FIAT/functional.py`), never hard-code indices.
+2. **Completion from functional type**: a normal-derivative node's completion partner is
+   the corresponding tangential-derivative node; the completion is a per-entity statement,
+   determined by which components of the derivative jet the dual basis lacks.


Ah, this can be done numerically -- you have, say, one direction and need to get a basis for the orthogonal complement...

rckirby · 2026-07-13T18:48:59Z

+   (FTC, quintic endpoint rule, integration by parts against trace moments) are all
+   instances of one computation — express a completion functional applied to the
+   polynomial space in the basis of the element's own nodes, i.e. solve with the
+   generalized Vandermonde matrix. This is where GEM-based dual evaluation comes in.


Does this GEM Vandermonde have to be done in physical space? If so, that's where we have a loss in performance (and perhaps eventually accuracy) over the hand-coded theory.

pbrubeck · 2026-07-13T20:01:58Z

+def make_unisolvent_points(element, interior=False):
+    degree = element.degree()
+    ref_complex = element.get_reference_complex()
+    top = ref_complex.get_topology()
+    pts = []
+    if interior:
+        dim = ref_complex.get_spatial_dimension()
+        for entity in top[dim]:
+            pts.extend(ref_complex.make_points(dim, entity, degree+dim+1, variant="gll"))
+    else:
+        for dim in top:
+            for entity in top[dim]:
+                pts.extend(ref_complex.make_points(dim, entity, degree, variant="gll"))
+    return pts
+
+
+def check_zany_mapping(element, ref_to_phys, *args, **kwargs):


restore these functions to their original location in test_zany_mapping.py

pbrubeck · 2026-07-13T20:08:01Z

+    L = numpy.einsum("jcq,c->jq", T.reshape(T.shape[0], -1, len(points)),
+                     ndir.ravel())


use numpy.tensordot(..., axes=)

pbrubeck · 2026-07-13T20:08:17Z

+        processed.add(i)
+
+
+def _piola_point_rows(V: numpy.ndarray, group: dict, J: Node,


migrate to the class

pbrubeck · 2026-07-13T20:08:33Z

+        raise NotImplementedError("Cannot yet Piola-transform this node group.")
+
+
+def _piola_facet_rows(V: numpy.ndarray, group: dict,


migrate to the class

rckirby

Few more notes

rckirby · 2026-07-13T21:31:13Z

+        raise NotImplementedError(
+            f"{type(self).__name__} does not implement automatic basis transformation.")
+
+    def _invariant_dofs(self, group, dim, sd):


Then this needs to be an @abstractmethod

rckirby · 2026-07-13T21:31:26Z

+        raise NotImplementedError(
+            f"{type(self).__name__} does not implement automatic basis transformation.")
+
+    def _facet_dof_rows(self, V, group, fiat_element, entity, J, processed):


@abstractmethod

rckirby · 2026-07-13T21:31:42Z

+        raise NotImplementedError(
+            f"{type(self).__name__} does not implement automatic basis transformation.")
+
+    def _point_dof_rows(self, V, group, fiat_element, J, processed):


@abstractmethod

rckirby · 2026-07-13T21:32:46Z

+            raise NotImplementedError(
+                f"{type(self).__name__} expects an affine pullback, not {mappings}.")
+
+    def _invariant_dofs(self, group, dim, sd):


If you have an integral moment with no derivatives, do you need to consider rescaling to get things matching up? In which case you have a non-unit diagonal entry.

numpy.linalg.pinv silently returns a least-squares fit if the group's tangential profiles are not actually independent; solving the small, square Gram system B @ B.T instead fails loudly on such a bug since B has full row rank by unisolvence. Document why the per-point frame coordinate profile (not just a tangent direction) is needed to tell facet dofs apart.

Replace the normal/tangential frame decomposition (FacetFrame, generalized_cross) with a direct duality argument: a physical node's row is B_ij = l_i(psi_hat_j), evaluated against the reference nodal basis with no geometric frame, by applying the chain rule to the tabulation and contracting the (untouched) weighted direction tensor. B is not V -- V = B^-1 still has to be assembled -- but B is block lower triangular by topological dimension, so V comes from inverting each entity's small diagonal block against the already-known lower-dimensional rows, reusing physically_mapped.inverse. Restrict the contraction to dofs with a nonzero tabulated derivative, following Walkington's sparsity pattern, since most of a nodal basis vanishes at any given handful of points. PiolaPhysicallyMappedElement is untouched.

…k/zany-auto

pbrubeck and others added 3 commits July 13, 2026 10:34

pbrubeck commented Jul 13, 2026

View reviewed changes

WIP

a865300

pbrubeck force-pushed the pbrubeck/zany-auto branch from d6e54fa to a865300 Compare July 13, 2026 11:52

pbrubeck mentioned this pull request Jul 13, 2026

TensorProduct(PhysicallyMapped) #258

Draft

pbrubeck and others added 8 commits July 13, 2026 15:11

Remove accidentally committed vim swap files

516ae87

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Enforce pydocstyle; document it in the PR expectations

296fb68

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

pbrubeck commented Jul 13, 2026

View reviewed changes

Comment thread finat/zany.py Outdated

pbrubeck commented Jul 13, 2026

View reviewed changes

Comment thread finat/zany.py

rckirby reviewed Jul 13, 2026

View reviewed changes

pbrubeck commented Jul 13, 2026

View reviewed changes

rckirby requested changes Jul 13, 2026

View reviewed changes

pbrubeck added 5 commits July 13, 2026 23:41

split AGENTS.md

4fa258b

split AGENTS.md

2d156f0

merge conflict

bddfcf4

Merge remote-tracking branch 'origin/pbrubeck/zany-auto' into pbrubec…

dbc95c5

…k/zany-auto

		class Functional:
		"""Symbolic degree of freedom with a single derivative direction.


		Key facts the framework rests on:

		* FIAT normals are "UFC consistent": computed from the tangents by the same formula

		L = numpy.einsum("jcq,c->jq", T.reshape(T.shape[0], -1, len(points)),
		ndir.ravel())

		processed.add(i)


		def _piola_point_rows(V: numpy.ndarray, group: dict, J: Node,

		raise NotImplementedError("Cannot yet Piola-transform this node group.")


		def _piola_facet_rows(V: numpy.ndarray, group: dict,

Uh oh!

Conversation

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

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!

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

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

rckirby left a comment

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!

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!

Choose a reason for hiding this comment

Uh oh!

rckirby left a comment

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!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

pbrubeck commented Jul 13, 2026 •

edited

Loading

pbrubeck Jul 13, 2026 •

edited

Loading