diff --git a/AGENTS.md b/AGENTS.md index af595278..dd2fa343 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -9,11 +9,7 @@ This document outlines the guidelines and architectural context for AI agents as When assisting with contributions to FIAT and the Firedrake project, AI agents and their human counterparts must adhere to the following strict policies: * The use of AI tools must be explicitly declared alongside the specific tool used. -* A human developer must lead the Pull Request. -* The human contributor must understand every change made to the codebase. * Reviewer questions must be answered directly by the human, rather than acting as a relay to the AI. -* Any generated code must be executed locally to verify that it functions correctly. -* AI tools must not be used to resolve issues that are labeled as 'good first issue'. --- @@ -28,6 +24,19 @@ FIAT operates within Firedrake's automated system for solving partial differenti --- +## Repository Layout + +Despite being called "fiat", this repository (Python package `firedrake_fiat`) contains **three** packages plus their tests, all in one tree: + +* `FIAT/` — reference-element definitions: cells, polynomial sets, dual bases (`FIAT/functional.py` holds the taxonomy of degree-of-freedom functionals), and element families. +* `finat/` — symbolic layer on top of FIAT. Physically mapped elements (`finat/physically_mapped.py`, `finat/argyris.py`, `finat/morley.py`, `finat/hct.py`, `finat/piola_mapped.py`, …) construct GEM expressions for basis transformations. See `finat/AGENTS.md` for package-specific implementation conventions. +* `gem/` — the tensor-algebra intermediate language used to express transformations symbolically (`gem/gem.py` for nodes, `gem/interpreter.py` to evaluate expressions numerically in tests). See `gem/AGENTS.md` for package-specific implementation conventions. +* `test/` — note the singular name: tests live at `test/finat/test_zany_mapping.py`, `test/FIAT/...`, `test/gem/...` (not `tests/`). + +Read the subdirectory `AGENTS.md` for any package you are about to edit — they hold +implementation-level conventions (GEM operator overloading rules, FInAT Jacobian sign +conventions, etc.) that this top-level file does not repeat. + ## Environment and Setup Bugs can exist within Firedrake or any of its component packages, explicitly including FIAT. To effectively develop and debug FIAT: @@ -52,6 +61,33 @@ Agents modifying FIAT code must follow these fundamental development principles: --- +## Style and Conventions + +When writing Python code for FIAT, maintain the ecosystem's structural and stylistic integrity: + +* Class attributes must be declared in one visible location. +* Attributes must be initialized in the `__init__` constructor or declared as a `functools.cached_property` if they are expensive to compute. +* Ad hoc lazy initialization discovering attributes via `hasattr`, `setattr`, or `getattr` scattered across methods is strictly prohibited. +* Boolean attributes must be used to record initialization intent and state instead of probing for the presence of state built by an initialization function. +* New code must include type hints on all function and method signatures. +* Public-facing APIs must include properly formatted `numpydoc`-style docstrings. + +--- + +## Pull Request Expectations + +* All changes are expected to arrive through GitHub Pull Requests. +* Keep diffs reviewable and focused. +* Before concluding work verify that the relevant subset of the pytest test suite + succeeds locally. +* CI (`.github/workflows/test.yml`) checks `flake8` and `pydocstyle .` (both + configured in `setup.cfg`); run them locally and fix all findings before pushing. +* Watch for `pydocstyle` D413 (blank line after the last numpydoc section) and D417 + (every argument, including keyword arguments, needs a description) in new docstrings. + clean `flake8` pass does not imply a clean `pydocstyle` pass. + +--- + ## Pattern Matching and Mathematical Reasoning When designing or debugging FIAT, FInAT, and GEM changes, use the existing codebase as a library of @@ -84,16 +120,76 @@ mathematical patterns rather than starting from ad hoc special cases: restriction, and tensor-product algebra agree with the element's continuity and approximation properties, not when one or two test cases happen to pass. -## Style and Conventions -When writing Python code for FIAT, maintain the ecosystem's structural and stylistic integrity: +## Pattern Matching and Mathematical Reasoning -* Class attributes must be declared in one visible location. -* Attributes must be initialized in the `__init__` constructor or declared as a `functools.cached_property` if they are expensive to compute. -* Ad hoc lazy initialization discovering attributes via `hasattr`, `setattr`, or `getattr` scattered across methods is strictly prohibited. -* Boolean attributes must be used to record initialization intent and state instead of probing for the presence of state built by an initialization function. -* New code must include type hints on all function and method signatures. -* Public-facing APIs must include properly formatted `numpydoc`-style docstrings. -* CI enforces `pydocstyle` (see the `[pydocstyle]` section of `setup.cfg` for the active ignore list) - in addition to `flake8`; run `pydocstyle ` locally before finishing a change, since a - clean `flake8` pass does not imply a clean `pydocstyle` pass. +This section used to be a rough draft of aspirations. It is now grounded in the concrete +experience of automating the Kirby (2017) / Aznaran-Kirby-Farrell (2022) / Brubeck & Kirby +(2025) transformation theory (`finat/zany.py`, `finat/functional.py`) for Morley, Hermite, +Argyris, Bell, Mardal-Tai-Winther, Johnson-Mercier, and Guzman-Neilan. The lessons below are +about *how to design and debug this kind of code*, not just about this one project. + +### Design strategies that generalize + +* **Design by duality, never by direct construction.** Do not try to write physical basis + functions directly (that requires inverting a Vandermonde system). Write the physical *node* + (a linear functional — point evaluation or moment, always computable), push it forward, and + solve for its expansion in the reference nodes. This is the master strategy behind the whole + framework and generalizes to any element whose nodes, not whose basis functions, are the + simple objects. +* **Recover mathematical type from data, not from a class hierarchy.** `from_fiat` never asks + "is this an `IntegralMomentOfNormalDerivative`?"; it reads `pt_dict`/`deriv_dict` and derives + (order, direction, rank) numerically. This makes the framework forward-compatible with FIAT + functional classes that do not exist yet, as long as they reduce to the same numeric shape — + which essentially all "linear functional built from point/derivative data" dofs do. Prefer + this kind of structural recovery over adding another `isinstance`/type-tag branch whenever a + new functional needs support. +* **Turn a mathematical case-split into a template method, not a conditional.** When a theorem + reads "for pullback type X do A, for pullback type Y do B" against an otherwise identical + algorithmic skeleton, encode the skeleton once (`ZanyPhysicallyMappedElement + .basis_transformation`) and the case-split as small, named hook methods on sibling mixins + (`ScalarPhysicallyMappedElement`, `PiolaPhysicallyMappedElement`). Keep the underlying linear + algebra as free functions taking plain arrays/GEM expressions with no `self` + (`FacetFrame`, `_scalar_facet_rows`, `_piola_facet_rows`, ...): the class layer should only + ever be a thin adapter from "element structure" to "which pure-math function to call," never + a place where new mathematics is derived. This also keeps the generic infrastructure + (`PhysicallyMappedElement` in `finat/physically_mapped.py`, used by hand-coded elements like + Arnold-Winther and HCT) free of any awareness of the theory-specific case-split. +* **Validate one element at a time against an independently computed ground truth, not just + against the closed-form answer you derived.** `test/finat/test_zany_mapping.py + ::check_zany_mapping` tabulates the FIAT element on an actual physical cell and least-squares + fits the physical values against the pulled-back reference tabulation — a computation + completely independent of the symbolic derivation. Its assertion message pretty-prints the + (row, column) location of the relative error; read those indices first, they usually + localize the bug to one entity/dof-type combination before you write any new code. +* **Generalize from one worked example, verify to machine precision, then extend — never + generalize speculatively ahead of a second concrete example.** The actual sequence that + worked: derive the mechanism on Morley alone, check it reproduces the deleted hand-coded + matrix bit-for-bit, *then* extend to Hermite, then Argyris/Bell, then the Piola family. Each + extension needed a genuinely new piece (interior invariance for Piola interior moments, + mixed-order jets for Bell/Argyris vertices, the reciprocal-basis correction for 3D) that a + premature generalization from Morley alone would not have anticipated. + +### Working with the human collaborator + +* **The literature is necessary but not sufficient — treat a stuck derivation as a cue to ask, + not to keep re-deriving.** Kirby (2017) and Aznaran-Kirby-Farrell (2022) work out 2D (or a + simplification that hides a 3D-only effect); none of the papers mention the reciprocal-basis + transformation law 3D Piola elements need. That fix came directly from the user's domain + expertise ("the papers do not consider the reciprocal basis that FIAT implements... Consider + this"), not from re-reading the papers harder. When a derivation matches every 2D case but + fails exactly where a paper's simplifying assumption would hide something, that is the moment + to ask what convention the *implementation* (not the paper) actually uses. +* **Expect, and design for, architectural correction mid-project.** This project's real + trajectory was: prototype on Morley → generalize into a `Functional`/free-function design → + "don't extend further, make the computation generic first, get rid of the very specific + helpers" → extend to Argyris/Bell → extend to the Piola family → "refactor into + `Scalar`/`PiolaPhysicallyMappedElement` mixins so the `if piola` disappears from the loop" → + "move the method body into the Mixin" (keep the generic `PhysicallyMappedElement` free of + zany-specific knowledge). None of these were bug fixes; each was the human reshaping the + design once its true shape became visible. Treat working code as a draft the architecture + will still move under, and re-run the full verification (test suite, `flake8`, `pydocstyle`) + after every such reshaping, not only after functional changes. +* **Keep this file's *why*, not just the *what*.** Git history has the diffs; this file exists + so the next session does not have to rediscover the reciprocal-basis subtlety, the + sign-invariance requirement, or the identity-row idiom by re-reading a diff from scratch. diff --git a/finat/AGENTS.md b/finat/AGENTS.md new file mode 100644 index 00000000..9ed7be9d --- /dev/null +++ b/finat/AGENTS.md @@ -0,0 +1,101 @@ +## FInAT implementation conventions + +* `coordinate_mapping.jacobian_at(point)` returns $\partial x_{\text{phys}} / \partial + x_{\text{ref}}$ — the **inverse** of the papers' $J$ (papers map physical → reference). + So where the paper writes $J^{-T}$, the code uses the FInAT Jacobian transposed; e.g. + Hermite's vertex blocks are `J[j, k]` entries of $M$ directly. +* `basis_transformation` builds the numpy object array `V` (row = reference node, column + = physical node, entries are GEM scalars) and returns `gem.ListTensor(V.T)` = $M$. + Start from `identity(ndof)` (in `finat/physically_mapped.py`) and overwrite the + non-identity rows. +* Existing generalized helpers, all in `finat/argyris.py`: + * `_jet_transform(J, order)`: chain-rule block for a symmetric derivative jet of any + order (order 1 → Jacobian, order 2 → $\Theta$, …), handling symmetric-component + flattening. + * `_vertex_transform(V, vorder, cell, mapping)`: places jet blocks for every vertex. + * `_normal_tangential_transform(cell, J, detJ, edge)`: returns $(B_{nn}, B_{nt}, Jt)$ + for an edge, expressing $B$ entries via $G^{TT}$ Gram data so only GEM-representable + quantities appear (`detJ/beta`, `alpha/beta`); 3D variant for faces is + `morley_transform` in `finat/walkington.py` (Morley itself is automated now; this + 3D helper survives only as Walkington's dependency). + * `_edge_transform(V, vorder, eorder, cell, mapping, avg)`: the Jacobi-moment edge rows + for integral-variant Argyris/HCT, encoding the endpoint values + $P_i^{(1,1)}(\pm 1)$ and the trace-moment coupling. +* Reduced/constrained elements (`ReducedHsiehCloughTocher` in `finat/hct.py`, `Bell`): + the FIAT element is the *extended* element (12 basis functions for reduced HCT), the + FInAT element exposes only $\nu$ dofs: `V = identity(numbf, ndof)` is rectangular, + `entity_dofs()` is overridden to empty the constrained entities, and + `space_dimension()` returns the reduced count. +* Conditioning convention: after assembling `V`, columns associated with derivative dofs + are rescaled by powers of `coordinate_mapping.cell_size()` (a per-vertex $h$). This + redefines the physical dofs as $h$-scaled derivatives — consistent across cells because + the scaling depends only on shared vertices/edges. Any new transformation must follow + the same convention or mass-matrix conditioning degrades + (`test/finat/test_mass_conditioning.py`). +* Verification: `test/finat/test_zany_mapping.py::check_zany_mapping` computes the exact + $M$ numerically by tabulating the FIAT element on a *physical* cell and least-squares + solving against pulled-back (and Piola-mapped, if applicable) reference tabulations, + then compares against `basis_transformation` evaluated through `gem.interpreter`. Its + assertion message pretty-prints the relative-error matrix with row/column indices — + read those indices to identify which dof couplings are wrong. + +## Mathematical structures to recognize in FIAT/FInAT + +* **A degree of freedom is fully described by five numbers, not by its FIAT class.** Every + functional FIAT builds from `pt_dict`/`deriv_dict` reduces to (points, weights, derivative + order $m$, a direction tensor of rank $m$, a value rank for vector/tensor-valued dofs). + `IntegralMomentOfNormalDerivative`, `PointNormalDerivative`, `TensorBidirectionalIntegralMoment`, + etc. are just different ways of *constructing* the same five numbers. Recognizing this + collapses "N functional types to support" into "one shape to recover numerically" + (`finat.PhysicallyMappedFunctional.from_fiat`, `finat/functional.py`). +* **Pullback is always "contract each tensor slot of the direction with a fixed matrix."** + Order-0 (values) are invariant (zero slots to contract). Order-$m$ derivatives contract $m$ + slots with the Jacobian $J$ (the chain rule). Rank-$r$ Piola values contract $r$ slots with + the cofactor matrix $K = \operatorname{adj}(J)^T$. This is *the same operation* with a + different matrix, which is why the scalar and Piola code in `finat/zany.py` are mirror images + of each other (`_scalar_point_rows` / `_piola_point_rows`, `_scalar_facet_rows` / + `_piola_facet_rows`) rather than unrelated implementations. +* **Frame decomposition is the one computational primitive underlying every non-affine + element.** Facet dofs (Morley/Argyris/Bell normal derivatives; MTW/JM/GN normal-tangential + moments) and vertex-jet completions (Hermite/Argyris/Bell gradients and Hessians) are all + solved the *same* way: split a direction/profile into an invariant part and a part that needs + completing, express the pulled-back quantity in the frame built from the mapped generators of + that split (`FacetFrame`, or the direction-basis inverse in `_scalar_point_rows`), and solve + symbolically via `adjugate`/`determinant`. Once this pattern is visible, "add a new element" + stops being "derive new math" and becomes "which invariant subspace, which frame." +* **Constrained/extended elements are the same construction as a restriction.** Bell and + Guzman-Neilan are both "take the extended FIAT element with its constraint functionals as + extra dofs, transform the whole thing, then keep only the first $\nu$ columns." This is + Kirby (2017) §5's extended-element proposition, and in code it is nothing more than + `space_dimension()` returning a smaller count than the FIAT element's — no special-casing + needed in the transformation loop itself. + +### Confusing mathematical ideas, clarified + +* **The Jacobian direction is flipped between the papers and the code.** The papers define + $F$ from physical to reference space, so their $J$ is FInAT's Jacobian *inverse*. + `coordinate_mapping.jacobian_at(point)` returns $\partial x_{\text{phys}}/\partial + x_{\text{ref}}$ (see "FInAT implementation conventions" below); a paper's $J^{-T}$ is FInAT's + plain Jacobian, transposed. This is a constant source of "why does my formula look + transposed" confusion — check this convention flip before suspecting a sign or index bug. +* **Components against a reciprocal/dual basis transform contragrediently — this is the single + subtlety that broke 3D and is not in any of the papers.** FIAT builds the tangential + component of 3D Piola-mapped dofs (MTW) using `cross(n, t_k)`, the *reciprocal* partner of + the tangent frame, not the tangent frame itself. A general fact from differential geometry: + if a frame transforms by a matrix $A$, components expressed against its *reciprocal* frame + transform by $A^{-T}$ (up to a determinant factor), not by $A$. Kirby (2017)/Aznaran-Kirby- + Farrell (2022) only work out the 2D case, where the tangent "plane" is 1-dimensional and this + correction $S^{-1}$ collapses to $1$ — silently hiding the effect. **Whenever a FIAT dual set + builds a direction via a cross product against the normal (a common way to get "the other + in-plane directions" in 3D), check whether its transformation law is contragredient before + assuming it matches the direct frame.** +* **Numerically recovered invariants (SVD, pinv) are only defined up to a group action, and + formulas built from them must be invariant under that action.** `from_fiat`'s SVD recovers a + direction/weight pair $(\hat n, w)$ up to a joint sign flip ($\hat n \to -\hat n$, $w \to + -w$ leaves the functional unchanged). A formula like $r_k = a x_k + (1-c)\beta_k$ that + implicitly assumes a particular sign will be wrong on exactly the subset of entities where + SVD happened to pick the other sign — a bug that looks like it "mostly works" and is + otherwise very hard to localize. The fix, $r_k = x_k - c\beta_k$, is invariant under the + joint flip. **Always test the full topology of a non-degenerate, non-symmetric physical + cell** (see `test/finat/conftest.py::MyMapping`'s deliberately irregular vertex coordinates) — + a symmetric test cell can accidentally hide a sign bug that only shows up on a generic mesh. diff --git a/finat/__init__.py b/finat/__init__.py index bde435a1..bcefc5f6 100644 --- a/finat/__init__.py +++ b/finat/__init__.py @@ -43,6 +43,7 @@ from .nodal_enriched import NodalEnrichedElement # noqa: F401 from .quadrature_element import QuadratureElement, make_quadrature_element # noqa: F401 from .restricted import RestrictedElement # noqa: F401 +from .functional import PhysicallyMappedFunctional # noqa: F401 from .runtime_tabulated import RuntimeTabulated # noqa: F401 from . import quadrature # noqa: F401 from . import cell_tools # noqa: F401 diff --git a/finat/alfeld_sorokina.py b/finat/alfeld_sorokina.py index c0ccf691..7eb7336f 100644 --- a/finat/alfeld_sorokina.py +++ b/finat/alfeld_sorokina.py @@ -13,7 +13,7 @@ def __init__(self, cell, degree=2): cite("AlfeldSorokina2016") super().__init__(FIAT.AlfeldSorokina(cell, degree)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): sd = self.cell.get_spatial_dimension() bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) diff --git a/finat/argyris.py b/finat/argyris.py index e907eff8..02f02766 100644 --- a/finat/argyris.py +++ b/finat/argyris.py @@ -139,7 +139,7 @@ def __init__(self, cell, degree=5, variant=None, avg=False): self.avg = avg super().__init__(fiat_element) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): sd = self.cell.get_spatial_dimension() top = self.cell.get_topology() diff --git a/finat/aw.py b/finat/aw.py index c9400031..df990e58 100644 --- a/finat/aw.py +++ b/finat/aw.py @@ -55,7 +55,7 @@ def __init__(self, cell, degree=2): cite("Arnold2003") super().__init__(FIAT.ArnoldWintherNC(cell, degree)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): """Note, the extra 3 dofs which are removed here correspond to the constraints.""" numbf = self._element.space_dimension() @@ -86,7 +86,7 @@ def __init__(self, cell, degree=3): cite("Arnold2002") super().__init__(FIAT.ArnoldWinther(cell, degree)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): # The extra 6 dofs removed here correspond to the constraints numbf = self._element.space_dimension() ndof = self.space_dimension() diff --git a/finat/bell.py b/finat/bell.py index d174fe16..aca22676 100644 --- a/finat/bell.py +++ b/finat/bell.py @@ -20,7 +20,7 @@ def __init__(self, cell, degree=5): reduced_dofs[sd-1][entity] = [] self._entity_dofs = reduced_dofs - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): # Jacobian at barycenter sd = self.cell.get_spatial_dimension() top = self.cell.get_topology() diff --git a/finat/functional.py b/finat/functional.py new file mode 100644 index 00000000..2fd75566 --- /dev/null +++ b/finat/functional.py @@ -0,0 +1,276 @@ +r"""Symbolic representation of degrees of freedom. + +A :class:`PhysicallyMappedFunctional` represents a degree of freedom in the form + +.. math:: \\ell(f) = \\sum_q w_q \\langle D, \\nabla^m f(x_q) \\rangle, + +where the points :math:`x_q` and quadrature/moment weights :math:`w_q` +are numeric, and the direction tensor :math:`D` (of rank equal to the +derivative order :math:`m`) may be numeric, for functionals defined on +the reference cell, or a GEM expression, for functionals carrying +physical geometry. + +This representation is the foundation for automating the transformation +theory of Kirby (2017): degrees of freedom of any FIAT element are +converted to this common form directly from their point and derivative +dictionaries, with the derivative direction recovered numerically, so +that no dispatch over FIAT functional types is required. + +The physical counterpart of a reference functional is assumed to share +its points and weights: integral moments must be measure-intrinsic +(e.g. integral averages), following the reference node convention of +Brubeck & Kirby (2025). +""" + +from functools import reduce +from math import factorial, prod +from operator import mul + +import numpy + +from FIAT.finite_element import FiniteElement +from FIAT.polynomial_set import mis +from FIAT.functional import Functional as FIATFunctional +from gem import Node, Zero + + +def multiindices(sd: int, order: int) -> list: + """Multi-indices of a given order, with axis ordering for order 1.""" + return sorted(mis(sd, order), reverse=True) + + +class PhysicallyMappedFunctional: + """Symbolic degree of freedom with a single derivative direction. + + Parameters + ---------- + points : + Tuple of reference-cell points. + weights : + Numeric weight for each point. + order : + The derivative order :math:`m`. + direction : + For ``order > 0``, the direction tensor of rank ``order``, + either numeric or a GEM expression; ``None`` for ``order == 0``. + divergence : + Whether this is a divergence functional: the trace, at each + point, of the tensor pairing an order-1 derivative multi-index + with a rank-1 value component. Unlike an ordinary derivative, + this does not have a single ``direction``, since the trace + pattern spans the full derivative and value index ranges; it is + recovered and handled as its own case because it commutes with + the (contravariant) Piola pullback up to the Jacobian + determinant, independently of the entity it sits on. + mapping : + The FIAT mapping string of the basis functions this functional + is dual to (``"affine"``, ``"contravariant piola"``, ...): the + type tag of the value slots, deciding which matrix each value + slot of the weight tensor is contracted with under push-forward, + exactly as the derivative slots of ``direction`` are contracted + with the Jacobian. + + """ + + def __init__(self, points: tuple, weights: numpy.ndarray, + order: int = 0, direction=None, rank: int = 0, + divergence: bool = False, mapping: str = "affine"): + self.points = points + self.weights = weights + self.order = order + self.direction = direction + self.rank = rank + self.divergence = divergence + self.mapping = mapping + + @classmethod + def from_fiat(cls, node: FIATFunctional, tol: float = 1e-12, + mapping: str = "affine") -> "PhysicallyMappedFunctional": + """Construct a symbolic PhysicallyMappedFunctional from a FIAT functional. + + The construction only inspects the point and derivative + dictionaries: the derivative order and the (common) direction of + differentiation are recovered numerically by factorizing the + matrix of derivative weights. + + Parameters + ---------- + node : + The FIAT functional. + tol : + Relative tolerance for the rank-one factorization of the + derivative weights. + mapping : + The FIAT mapping string of the basis functions this + functional is dual to. + + Returns + ------- + PhysicallyMappedFunctional + The symbolic representation of the FIAT functional. + + """ + if node.pt_dict and node.deriv_dict: + raise NotImplementedError( + f"{type(node).__name__} mixes value and derivative weights.") + + if not node.deriv_dict: + points = tuple(node.pt_dict) + comps = {comp for pt in points for w, comp in node.pt_dict[pt]} + rank = len(max(comps)) + if rank == 0: + weights = numpy.asarray([w for pt in points + for w, comp in node.pt_dict[pt]]) + return cls(points, weights, mapping=mapping) + # value weight profile: one row of component weights per point + sd = node.ref_el.get_spatial_dimension() + weights = numpy.zeros((len(points), sd**rank)) + shape = (sd,) * rank + for q, pt in enumerate(points): + for w, comp in node.pt_dict[pt]: + weights[q, numpy.ravel_multi_index(comp, shape)] += w + return cls(points, weights, rank=rank, mapping=mapping) + + sd = node.ref_el.get_spatial_dimension() + order = node.max_deriv_order + alphas = multiindices(sd, order) + lookup = {alpha: k for k, alpha in enumerate(alphas)} + + points = tuple(node.deriv_dict) + has_comp = any(comp != tuple() for pt in points + for w, alpha, comp in node.deriv_dict[pt]) + if has_comp: + # A divergence: at each point, the weights pairing an + # order-1 derivative multi-index with a rank-1 value + # component must form a scalar multiple of the trace. + if order != 1: + raise NotImplementedError( + f"{type(node).__name__} has vector components.") + weights = numpy.zeros(len(points)) + for q, pt in enumerate(points): + Wq = numpy.zeros((sd, sd)) + for w, alpha, comp in node.deriv_dict[pt]: + if len(comp) != 1: + raise NotImplementedError( + f"{type(node).__name__} has vector components.") + Wq[lookup[tuple(alpha)], comp[0]] += w + if not numpy.allclose(Wq, Wq[0, 0] * numpy.eye(sd), atol=tol): + raise NotImplementedError( + f"{type(node).__name__} is not a divergence functional.") + weights[q] = Wq[0, 0] + return cls(points, weights, order=order, divergence=True, + mapping=mapping) + + W = numpy.zeros((len(points), len(alphas))) + for q, pt in enumerate(points): + for w, alpha, comp in node.deriv_dict[pt]: + W[q, lookup[tuple(alpha)]] += w + + # Factor the weights as a common direction times scalar weights + u, s, vt = numpy.linalg.svd(W) + if any(s[1:] > tol * s[0]): + raise NotImplementedError( + f"{type(node).__name__} has no common derivative direction.") + direction = vt[0] + weights = u[:, 0] * s[0] + return cls(points, weights, order=order, direction=direction, + mapping=mapping) + + def with_direction(self, direction) -> "PhysicallyMappedFunctional": + """Return the same functional with another direction tensor.""" + return type(self)(self.points, self.weights, + order=self.order, direction=direction, + mapping=self.mapping) + + def pullback(self, A) -> "PhysicallyMappedFunctional": + r"""Contract each derivative slot of the direction with a matrix. + + By the chain rule, reference derivatives of a pullback are + physical derivatives contracted with the Jacobian, so viewing + this reference functional as acting on physical functions + amounts to contracting each slot of the direction tensor with + the cell Jacobian; more general per-slot matrices arise when + the direction is instead mapped through a physical frame. + + Parameters + ---------- + A : + The per-slot matrix: a GEM expression of shape ``(sd, sd)``, + or an ``(sd, sd)`` array with numeric or GEM scalar entries. + + Returns + ------- + PhysicallyMappedFunctional + The functional with direction :math:`A \\otimes \\dots + \\otimes A : D`, collapsed back onto multi-index + coefficients. + + """ + if self.order == 0: + return self + if isinstance(A, Node): + sd = A.shape[0] + A = numpy.array([[A[i, k] for k in range(sd)] for i in range(sd)], + dtype=object) + A = numpy.asarray(A) + sd = A.shape[0] + symbolic = A.dtype == object + alphas = multiindices(sd, self.order) + lookup = {alpha: k for k, alpha in enumerate(alphas)} + + direction = numpy.full(len(alphas), Zero() if symbolic else 0.0, + dtype=object if symbolic else float) + shape = (sd,) * self.order + for index in numpy.ndindex(shape): + # Distribute the multi-index coefficients over a symmetric tensor + alpha = _index_alpha(index, sd) + scale = prod(map(factorial, alpha)) / factorial(self.order) + d = self.direction[lookup[alpha]] * scale + if not isinstance(d, Node) and d == 0: + continue + # Contract each slot with A, collapsing onto multi-indices + for target in numpy.ndindex(shape): + term = reduce(mul, (A[target[k], index[k]] + for k in range(self.order)), d) + m = lookup[_index_alpha(target, sd)] + direction[m] = direction[m] + term + return self.with_direction(direction) + + def evaluate(self, fiat_element: FiniteElement) -> numpy.ndarray: + r"""Apply this functional to the nodal basis of a FIAT element. + + This is the generalized Vandermonde computation: the restriction + of a functional :math:`\\ell` to the polynomial space satisfies + :math:`\\pi \\ell = \\sum_j \\ell(\\psi_j)\\, \\pi n_j`. Only + valid for functionals with numeric direction. + + Parameters + ---------- + fiat_element : + The FIAT element providing the nodal basis. + + Returns + ------- + numpy.ndarray + The vector of values of this functional on the nodal basis. + + """ + sd = fiat_element.get_reference_element().get_spatial_dimension() + tab = fiat_element.tabulate(self.order, self.points) + if self.rank > 0: + T = tab[(0,) * sd] + T = T.reshape(T.shape[0], -1, len(self.points)) + return numpy.einsum("jcq,qc->j", T, self.weights) + if self.order == 0: + return tab[(0,) * sd] @ self.weights + alphas = multiindices(sd, self.order) + return self.direction @ numpy.array([tab[alpha] @ self.weights + for alpha in alphas]) + + +def _index_alpha(index: tuple, sd: int) -> tuple: + """Convert a tensor index into a derivative multi-index.""" + alpha = [0] * sd + for k in index: + alpha[k] += 1 + return tuple(alpha) diff --git a/finat/hermite.py b/finat/hermite.py index 5b3a2811..fd5f6d76 100644 --- a/finat/hermite.py +++ b/finat/hermite.py @@ -11,7 +11,7 @@ def __init__(self, cell, degree=3): cite("Ciarlet1972") super().__init__(FIAT.CubicHermite(cell)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): Js = [coordinate_mapping.jacobian_at(vertex) for vertex in self.cell.get_vertices()] diff --git a/finat/hz.py b/finat/hz.py index 4e84da32..33e751dd 100644 --- a/finat/hz.py +++ b/finat/hz.py @@ -13,7 +13,7 @@ def __init__(self, cell, degree=3, variant=None, quad_scheme=None): self.variant = variant super().__init__(FIAT.HuZhang(cell, degree, variant=variant, quad_scheme=quad_scheme)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): ndofs = self.space_dimension() V = identity(ndofs) diff --git a/finat/johnson_mercier.py b/finat/johnson_mercier.py index 7b6d485c..7005cc93 100644 --- a/finat/johnson_mercier.py +++ b/finat/johnson_mercier.py @@ -13,7 +13,7 @@ def __init__(self, cell, degree=1, variant=None, quad_scheme=None): self._indices = slice(None, None) super().__init__(FIAT.JohnsonMercier(cell, degree, variant=variant, quad_scheme=quad_scheme)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): numbf = self._element.space_dimension() ndof = self.space_dimension() diff --git a/finat/morley.py b/finat/morley.py index ca3e89e8..fbd8655f 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -37,7 +37,7 @@ def __init__(self, cell, degree=2): cite("MingXu2006") super().__init__(FIAT.Morley(cell, degree=degree)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): sd = self.cell.get_spatial_dimension() top = self.cell.get_topology() # Jacobians at barycenter diff --git a/finat/mtw.py b/finat/mtw.py index 3e79c5d0..4086fa68 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -16,7 +16,7 @@ def __init__(self, cell, order=1): cite("Xie2008") super().__init__(FIAT.MardalTaiWinther(cell, order=order)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): sd = self.cell.get_spatial_dimension() bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py index 511d3e14..ecea01b0 100644 --- a/finat/physically_mapped.py +++ b/finat/physically_mapped.py @@ -1,10 +1,65 @@ +r"""Physically mapped elements and the automatic basis transformation. + +The transformation matrix of a physically mapped element is obtained by +duality alone (Kirby 2017, Brubeck & Kirby 2025). With +:math:`\hat\psi_j` the reference nodal basis, :math:`F` the cell map and +:math:`n_i` the physical node, + +.. math:: B_{ij} = n_i(\hat\psi_j \circ F^{-1}), \qquad V = B^{-1}, + +and the physical basis functions are :math:`M F^*(\hat\Psi)` with +:math:`M = V^T`. Two ingredients make the rows of :math:`B` computable +without frame algebra, for any element whose dual basis parses into +:class:`~finat.functional.PhysicallyMappedFunctional`: + +* **Physical nodes by per-slot maps.** The physical node shares the + points and weights of its reference partner; only the directional + data changes, one tensor slot at a time. A derivative slot of a + facet node maps its unit-normal component to the unit physical normal + :math:`K\hat{n}/|K\hat{n}|` (:math:`K = \operatorname{adj}(J)^T` the + cofactor matrix, which maps normals to normals) and its tangential + complement by :math:`J` (mapped tangents); away from facets + derivative slots keep their reference (Cartesian) directions. A + contravariant value slot of a facet moment maps its scaled-normal + component by :math:`K` (the cofactor lemma :math:`K\hat\nu^s = + \nu^s` is exact) and its tangential complement by the reciprocal of + the scalar tangential push-forward; Cartesian point data keeps its + weights, interior moments are invariant by convention, and + divergence nodes contract to :math:`\det J` times the identity. + +* **The adjoint acts on the tabulation.** Dually to the push-forward + law :math:`d = J^{\otimes m}\hat{d}`, each derivative slot of the + *numeric* reference tabulation carries :math:`J^{-T}` and each value + slot :math:`J/\det J`; transposing these onto the direction gives a + single effective per-slot matrix, so a row of :math:`B` is a numeric + generalized Vandermonde pairing with symbolic per-slot coefficients. + +Crucially, :math:`B` is never inverted densely, and its sparsity is +*inferred from numerical tabulations*: :math:`B` is assembled +numerically at a few generic sample Jacobians, entries that vanish at +every sample are structural zeros (the support law: dof :math:`j` +couples to dof :math:`i` only on the closure of :math:`i`'s entity, the +same property that makes the element conforming), and entries that +agree at every sample are constants, so push-forward-invariant rows +never touch symbolic algebra. The coupling graph is block +triangularized by its strongly connected components and :math:`V` is +computed by sparse block back-substitution, one small +adjugate/determinant solve per diagonal block, with fill-in confined to +the closures. +""" + from abc import ABCMeta, abstractmethod from collections.abc import Mapping +from functools import reduce +from operator import add, mul import gem import numpy +from FIAT.finite_element import FiniteElement + from finat.citations import cite +from finat.functional import PhysicallyMappedFunctional, multiindices class NeedsCoordinateMappingElement(metaclass=ABCMeta): @@ -65,7 +120,27 @@ def __len__(self): class PhysicallyMappedElement(NeedsCoordinateMappingElement): """A mixin that applies a "physical" transformation to tabulated - basis functions.""" + basis functions. + + :meth:`basis_transformation` derives the transformation generically + from the FIAT dual basis by the duality formulation of the module + docstring; the mixins of :mod:`finat.zany` supply the + family-specific conventions (expected pullback, moment + normalization, dof rescaling), and elements whose dual basis the + generic engine cannot parse override it entirely by hand. + """ + + #: Numerical tolerance used to drop negligible terms of the numeric + #: tabulations from the symbolic rows. + tol = 1e-12 + + #: Row-relative tolerance deciding, from the sample assemblies, which + #: entries of the transformation are structural zeros or constants. + pattern_tol = 1e-10 + + #: If False, physical scalar facet moments are plain integrals rather + #: than the measure-intrinsic integral averages of the reference nodes. + avg = True def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) @@ -73,13 +148,199 @@ def __init__(self, *args, **kwargs): cite("Kirby2019zany") self.restriction_indices = None - @abstractmethod - def basis_transformation(self, coordinate_mapping): - """Transformation matrix for the basis functions. + def basis_transformation(self, coordinate_mapping) -> gem.ListTensor: + r"""Transformation matrix for the basis functions. + + Assembles the generalized Vandermonde matrix :math:`B_{ij} = + n_i(\hat\psi_j\circ F^{-1})` and computes :math:`M = B^{-T}` + by sparse block back-substitution. The sparsity and the + constant entries of :math:`B` are inferred from numeric + assemblies at sample Jacobians; only the remaining entries are + built symbolically, in numerator/denominator form, and each + diagonal block of the strongly-connected-component + triangularization is solved through its adjugate. + + Parameters + ---------- + coordinate_mapping : + Object providing the physical geometry as GEM expressions. + + Returns + ------- + gem.ListTensor + The transformation matrix, shape ``(space_dimension(), + nbf)``, with the trailing constraint columns of an extended + element truncated. - :arg coordinate_mapping: Object providing physical geometry.""" + """ + fiat_element = self._element + self._check_mapping(fiat_element) + ref_el = fiat_element.get_reference_element() + sd = ref_el.get_spatial_dimension() + bary, = ref_el.make_points(sd, 0, sd + 1) + J = _materialize_jacobian(coordinate_mapping.jacobian_at(bary)) + + nodes = fiat_element.dual_basis() + mappings = fiat_element.mapping() + entity_ids = fiat_element.entity_dofs() + nbf = fiat_element.space_dimension() + ndof = self.space_dimension() + + # Parse the dual basis. FIAT may list a dof on more than one + # entity (e.g. the edge moments of Arnold-Winther reappear in its + # interior list); the lowest-dimensional entity owns the dof. + # Unparseable constraint functionals of an extended element are + # never exposed as physical dofs: their physical counterparts are + # defined as the pullbacks of the reference ones (Kirby 2017, + # section 5), keeping identity rows, and the corresponding + # columns are truncated below. + ells = {} + owners = {} + for dim in sorted(entity_ids): + for entity in sorted(entity_ids[dim]): + for i in entity_ids[dim][entity]: + if i in owners: + continue + owners[i] = (dim, entity) + try: + ells[i] = PhysicallyMappedFunctional.from_fiat( + nodes[i], mapping=mappings[i]) + except NotImplementedError: + if i < ndof: + raise + ells[i] = None + + # Infer the sparsity and the constant entries numerically + tabs = {} + samples = [self._assemble_numeric(fiat_element, ells, owners, Jk, tabs) + for Jk in _sample_jacobians(sd)] + pattern, constant, Bconst = _infer_pattern(samples, self.pattern_tol) + + # Build the non-constant rows symbolically, in + # numerator/denominator form, on their own pattern columns only + numer = numpy.empty(nbf, dtype=object) + denom = numpy.empty(nbf, dtype=object) + for i in range(nbf): + cols = numpy.flatnonzero(pattern[i]) + if ells[i] is None or constant[i, cols].all(): + numer[i] = numpy.array([_as_gem(v) for v in Bconst[i, cols]], + dtype=object) + denom[i] = one + else: + dim, entity = owners[i] + numer[i], denom[i] = _node_row( + fiat_element, ells[i], dim, entity, J, + self.avg, self.tol, tabs, cols=cols) + + # Sparse block back-substitution over the strongly connected + # components of the coupling graph, in dependency order: + # V_S = B_SS^{-1} (I_S - B_Sc V_c), with B = diag(1/denom) numer, + # so the denominators only enter the right-hand side. + V = numpy.full((nbf, nbf), zero, dtype=object) + inv_cache = {} + for block in _scc_blocks(pattern): + inside = set(block) + rhs = numpy.full((len(block), nbf), zero, dtype=object) + BSS = numpy.full((len(block), len(block)), zero, dtype=object) + index = {i: k for k, i in enumerate(block)} + for k, i in enumerate(block): + rhs[k, i] = denom[i] + for c, j in enumerate(numpy.flatnonzero(pattern[i])): + if j in inside: + BSS[k, index[j]] = numer[i][c] + elif not isinstance(numer[i][c], gem.Zero): + # numpy elementwise product: the object array must + # be the left operand, lest GEM broadcast a vector + rhs[k] = rhs[k] - V[j] * numer[i][c] + V[block] = _solve_block(BSS, rhs, inv_cache) + + self._rescale_dofs(V, fiat_element, coordinate_mapping) + return gem.ListTensor(V[:, :ndof].T) + + def _assemble_numeric(self, fiat_element: FiniteElement, ells: dict, + owners: dict, J: numpy.ndarray, tabs: dict) -> numpy.ndarray: + """Assemble B numerically at one sample Jacobian. + + :arg fiat_element: The FIAT element defined on the reference cell. + :arg ells: Map from dof index to parsed + :class:`~finat.functional.PhysicallyMappedFunctional`, or + None for a constraint functional keeping an identity row. + :arg owners: Map from dof index to its owning ``(dim, entity)``. + :arg J: Numeric sample Jacobian. + :arg tabs: Cache of reference tabulations, shared across samples. + :returns: The matrix B as a numpy array. + """ + B = numpy.eye(fiat_element.space_dimension()) + for i, ell in ells.items(): + if ell is None: + continue + dim, entity = owners[i] + res = _node_row(fiat_element, ell, dim, entity, J, + self.avg, self.tol, tabs) + if res is not None: + row, den = res + B[i] = row / den + return B + + def _check_mapping(self, fiat_element: FiniteElement) -> None: + """Verify that this class knows how to transform this element's pullback. + + The generic engine dispatches on the mapping of each dof; this + hook lets the family mixins of :mod:`finat.zany` reject FIAT + elements outside their conventions up front. + + :arg fiat_element: The FIAT element defined on the reference cell. + :raises NotImplementedError: If the pullback is not supported. + """ pass + def _rescale_dofs(self, V: numpy.ndarray, fiat_element: FiniteElement, + coordinate_mapping) -> None: + r"""Rescale the physical degrees of freedom by powers of the cell size. + + Each column of ``V`` is multiplied by :meth:`dof_scale` evaluated + with the cell size averaged over the vertices of the dof's entity. + This is the FInAT convention keeping the mass matrix + well-conditioned; it is consistent across cells because the + scaling only depends on shared entities. + + :arg V: Object array being assembled; columns are rescaled in place. + :arg fiat_element: The FIAT element defined on the reference cell. + :arg coordinate_mapping: Object providing the physical geometry as + GEM expressions. + """ + # cell_size may be a GEM expression or a numpy array of numbers + h = coordinate_mapping.cell_size() + top = fiat_element.get_reference_element().get_topology() + nodes = fiat_element.dual_basis() + entity_ids = fiat_element.entity_dofs() + for dim in entity_ids: + for entity in entity_ids[dim]: + verts = top[dim][entity] + havg = reduce(add, (h[v] for v in verts)) / len(verts) + for i in entity_ids[dim][entity]: + scale = self.dof_scale(nodes[i], dim, havg) + if scale is not None: + V[:, i] = V[:, i] * scale + + def dof_scale(self, node, dim: int, havg): + r"""Return the conditioning rescaling factor of one physical dof. + + The default convention redefines each physical node of derivative + order :math:`m > 0` with a factor :math:`h^{-m}`; elements whose + hand-coded transformations established a different convention + (e.g. the :math:`h^{-2}` vertex values of Hu-Zhang) override this + method. + + :arg node: The FIAT functional of the dof. + :arg dim: Topological dimension of the entity the dof sits on. + :arg havg: GEM scalar for the cell size averaged over the vertices + of the dof's entity. + :returns: The GEM scaling factor, or ``None`` for no rescaling. + """ + order = node.max_deriv_order + return havg**(-order) if order > 0 else None + def map_tabulation(self, ref_tabulation, coordinate_mapping): assert coordinate_mapping is not None M = self.basis_transformation(coordinate_mapping) @@ -277,3 +538,352 @@ def inverse(A): M[numpy.ix_(ids, ids)] = Minv return M + + +def _materialize_jacobian(J: gem.Node) -> numpy.ndarray: + """Materialize a GEM Jacobian node as a numpy object array of GEM scalars. + + :arg J: GEM expression for the cell Jacobian, shape ``(sd, sd)``. + :returns: A ``(sd, sd)`` numpy object array with the same entries as + ``J``, usable with the numeric-or-symbolic linear algebra of this + module (:func:`determinant`, :func:`adjugate`), which index a GEM + node's entries directly rather than through numpy fancy indexing. + """ + sd = J.shape[0] + return numpy.array([[J[i, k] for k in range(sd)] for i in range(sd)], dtype=object) + + +def _as_gem(v: float) -> gem.Node: + """Convert a number to a GEM constant, preserving structural zeros. + + :arg v: The number. + :returns: ``Zero`` for an exact zero, so that downstream products + and sums fold, the shared ``one`` for an exact one, and a + ``Literal`` otherwise. + """ + if v == 0: + return zero + if v == 1: + return one + return gem.Literal(v) + + +def _as_gem_array(A) -> numpy.ndarray: + """Convert a numeric array to an object array of GEM constants. + + :arg A: The numeric array. + :returns: An object array of the same shape, entrywise :func:`_as_gem`. + """ + A = numpy.asarray(A) + out = numpy.full(A.shape, zero, dtype=object) + for index in numpy.ndindex(A.shape): + out[index] = _as_gem(float(A[index])) + return out + + +def _sample_jacobians(sd: int, nsamples: int = 3) -> list: + """Deterministic generic sample Jacobians for sparsity inference. + + :arg sd: The spatial dimension. + :arg nsamples: How many samples; the last is orientation-reversing. + :returns: A list of well-conditioned ``(sd, sd)`` numeric Jacobians + with no special alignments, at which entries of the + transformation that vanish (or agree) at every sample may be + declared structurally zero (or constant). + """ + samples = [] + t = 0.0 + while len(samples) < nsamples: + t += 1.0 + A = numpy.eye(sd) + 0.4 * numpy.cos( + t + 7.3 * numpy.arange(1.0, sd * sd + 1.0)).reshape(sd, sd) + det = numpy.linalg.det(A) + if abs(det) < 0.3: + continue + if len(samples) == nsamples - 1 and det * numpy.linalg.det(samples[0]) > 0: + A = A[::-1] if sd > 1 else -A + samples.append(A) + return samples + + +def _infer_pattern(samples: list, tol: float) -> tuple: + """Infer the sparsity and the constant entries of B from samples. + + :arg samples: Numeric assemblies of B at the sample Jacobians. + :arg tol: Row-relative tolerance below which an entry is a + structural zero, and spread below which it is a constant. + :returns: A tuple ``(pattern, constant, Bconst)`` of the boolean + coupling pattern (with full diagonal), the boolean mask of + J-independent entries, and their values (averaged over the + samples, snapped to exact integers where they round). + """ + stack = numpy.stack(samples) + scale = numpy.abs(stack).max(axis=(0, 2)) + pattern = (numpy.abs(stack) > tol * scale[None, :, None]).any(axis=0) + numpy.fill_diagonal(pattern, True) + Bconst = stack.mean(axis=0) + spread = stack.max(axis=0) - stack.min(axis=0) + constant = spread <= tol * numpy.maximum(1.0, numpy.abs(Bconst)) + snapped = numpy.round(Bconst) + near = constant & (numpy.abs(Bconst - snapped) + <= tol * numpy.maximum(1.0, numpy.abs(Bconst))) + Bconst[near] = snapped[near] + return pattern, constant, Bconst + + +def _scc_blocks(pattern: numpy.ndarray) -> list: + """Strongly connected components of the coupling graph, in dependency order. + + Tarjan's algorithm on the digraph with an edge ``i -> j`` whenever + ``pattern[i, j]``: each emitted component only depends on (has + pattern entries in) previously emitted ones, so emission order is + elimination order, and the block triangular structure of B is a + consequence of the inferred sparsity rather than an assumption. + + :arg pattern: Boolean coupling pattern of B. + :returns: List of sorted lists of dof indices, one per component. + """ + n = pattern.shape[0] + adjacency = [numpy.flatnonzero(pattern[i]) for i in range(n)] + index = [-1] * n + low = [0] * n + onstack = [False] * n + stack = [] + blocks = [] + counter = [0] + + def strongconnect(v): + index[v] = low[v] = counter[0] + counter[0] += 1 + stack.append(v) + onstack[v] = True + for w in adjacency[v]: + if w == v: + continue + if index[w] == -1: + strongconnect(w) + low[v] = min(low[v], low[w]) + elif onstack[w]: + low[v] = min(low[v], index[w]) + if low[v] == index[v]: + component = [] + while True: + w = stack.pop() + onstack[w] = False + component.append(w) + if w == v: + break + blocks.append(sorted(component)) + + for v in range(n): + if index[v] == -1: + strongconnect(v) + return blocks + + +def _solve_block(B: numpy.ndarray, rhs: numpy.ndarray, cache: dict) -> numpy.ndarray: + """Solve one diagonal block through its adjugate. + + :arg B: Object array with the (GEM) block of B numerators. + :arg rhs: Object array with the right-hand side rows. + :arg cache: Cache of adjugate/determinant pairs, shared across + blocks so that identical blocks (e.g. the vertex jets of every + vertex) are only inverted once. + :returns: The rows of V for this block. + """ + if B.shape[0] == 1: + d = B[0, 0] + return rhs if d == one else rhs / d + key = gem.ListTensor(B) + try: + adj, det = cache[key] + except KeyError: + adj, det = cache.setdefault(key, (adjugate(B), determinant(B))) + return (adj @ rhs) / det + + +def _tabulate(fiat_element: FiniteElement, order: int, points: tuple, + tabs: dict) -> dict: + """Cached reference tabulation of the nodal basis. + + :arg fiat_element: The FIAT element. + :arg order: The derivative order. + :arg points: The points, as a hashable tuple. + :arg tabs: The cache. + :returns: The FIAT tabulation dict. + """ + key = (order, points) + try: + return tabs[key] + except KeyError: + return tabs.setdefault(key, fiat_element.tabulate(order, points)) + + +def _slot_map(ref_el, entity: int, J: numpy.ndarray, + derivative: bool, facet: bool) -> tuple: + r"""Effective per-slot matrix pairing a physical node with the reference tabulation. + + Transposing the per-slot action of the push-forward onto the + physical direction data gives a single matrix per slot: + :math:`G = J^{-1}\Phi` for a derivative slot with physical + direction map :math:`\Phi`, and :math:`G = (J/\det J)^T\Phi` for a + contravariant value slot with physical weight map :math:`\Phi`. + Away from a facet :math:`\Phi` is the identity (Cartesian data); + on a facet it maps the normal component by the cofactor law and the + tangential complement by the mapped tangents, as in the module + docstring. The map is returned in numerator/denominator form so + that the numerators stay free of symbolic division. + + :arg ref_el: The reference cell. + :arg entity: The facet number (ignored unless ``facet``). + :arg J: The cell Jacobian: numeric array, or object array of GEM scalars. + :arg derivative: Whether the slot is a derivative or a value slot. + :arg facet: Whether the node data is in the facet frame rather than + Cartesian. + :returns: The pair ``(numerator, denominator)`` with + ``G = numerator / denominator``. + """ + detJ = determinant(J) + adjJ = adjugate(J) + if not facet: + return (adjJ, detJ) if derivative else (J.T, detJ) + sd = ref_el.get_spatial_dimension() + symbolic = J.dtype == object + lit = _as_gem_array if symbolic else numpy.asarray + K = adjJ.T + if derivative: + n = ref_el.compute_normal(entity) + nn = n @ n + Kn = K @ n + q = (Kn @ Kn) ** 0.5 + Gnum = (numpy.outer(adjJ @ Kn, lit(n * (nn ** 0.5 / nn))) + + lit(numpy.eye(sd) - numpy.outer(n, n) / nn) * (q * detJ)) + return Gnum, detJ * q + nu = ref_el.compute_scaled_normal(entity) + nn = nu @ nu + Knu = K @ nu + qq = Knu @ Knu + IP = lit((numpy.eye(sd) - numpy.outer(nu, nu) / nn) / nn) + Gnum = (lit(numpy.outer(nu, nu) / nn) * (detJ * detJ) + + ((J.T @ K) * qq - numpy.outer(lit(nu), K.T @ Knu) * detJ) @ IP) + return Gnum, detJ * detJ + + +def _sparse_combination(coeffs, N: numpy.ndarray, tol: float, + symbolic: bool) -> numpy.ndarray: + """Combine numeric tabulation rows with (possibly symbolic) coefficients. + + Computes ``row[j] = sum_k coeffs[k] * N[k, j]``, dropping negligible + numeric entries so that whether a coupling is present is decided + from the numeric array, never by inspecting a symbolic expression. + + :arg coeffs: The coefficients, numeric or GEM scalars. + :arg N: The numeric tabulation pairing matrix. + :arg tol: Relative tolerance below which entries of ``N`` are dropped. + :arg symbolic: Whether to accumulate GEM expressions. + :returns: The combined row. + """ + row = (numpy.full(N.shape[1], zero, dtype=object) if symbolic + else numpy.zeros(N.shape[1])) + scale = numpy.abs(N).max() + if scale == 0: + return row + for k, c in enumerate(coeffs): + if isinstance(c, gem.Zero) or (not isinstance(c, gem.Node) and c == 0): + continue + for j in numpy.flatnonzero(numpy.abs(N[k]) > tol * scale): + row[j] = row[j] + c * N[k, j] + return row + + +def _node_row(fiat_element: FiniteElement, ell: PhysicallyMappedFunctional, + dim: int, entity: int, J: numpy.ndarray, avg: bool, tol: float, + tabs: dict, cols=None): + r"""One row of the generalized Vandermonde matrix B. + + The row is the pairing of the physical node data with the reference + tabulation carrying the adjoint of the push-forward, assembled with + the single effective per-slot matrix of :func:`_slot_map`; + divergence nodes commute with the Piola pullback up to + :math:`\det J`, and nodes whose physical counterpart is the + push-forward of the reference one by convention (point values, + interior moments) return None for an identity row. + + :arg fiat_element: The FIAT element. + :arg ell: The parsed reference node. + :arg dim: The dimension of the entity the node sits on. + :arg entity: The entity number. + :arg J: The cell Jacobian: numeric array, or object array of GEM scalars. + :arg avg: Whether physical scalar facet moments are integral averages. + :arg tol: Relative tolerance below which numeric couplings are dropped. + :arg cols: Optional column indices to restrict the row to. + :returns: The pair ``(numerator_row, denominator)`` with + ``B[i] = numerator_row / denominator``, or None for a row of + the identity. + """ + ref_el = fiat_element.get_reference_element() + sd = ref_el.get_spatial_dimension() + symbolic = J.dtype == object + if cols is None: + cols = slice(None) + + if ell.divergence: + tab = _tabulate(fiat_element, 1, ell.points, tabs) + alphas = [tuple(int(k == c) for k in range(sd)) for c in range(sd)] + div = sum(tab[alphas[c]].reshape(tab[alphas[c]].shape[0], sd, -1)[:, c, :] + for c in range(sd)) + row = div @ ell.weights + scale = numpy.abs(row).max() + row[numpy.abs(row) <= tol * scale] = 0.0 + row = row[cols] + if symbolic: + row = numpy.array([_as_gem(v) for v in row], dtype=object) + return row, determinant(J) + + if ell.rank: + if ell.mapping not in ("contravariant piola", "double contravariant piola"): + raise NotImplementedError( + f"Cannot transform value nodes with a {ell.mapping} pullback.") + if ell.order: + raise NotImplementedError( + "Cannot transform nodes mixing values and derivatives.") + npts = len(ell.points) + # Single-point rank-1 nodes are Cartesian point data, not facet + # moments, even when they sit on a codimension-1 entity (e.g. the + # edge midpoint values of Alfeld-Sorokina); single-point rank-2 + # nodes remain frame dofs (e.g. the edge dofs of the Hu-Zhang + # point variant). + point_data = npts == 1 and (ell.rank == 1 or dim != sd - 1) + if dim == sd and not point_data: + # Interior moments are Piola invariant by convention: their + # physical test functions are themselves Piola-mapped. + return None + facet = dim == sd - 1 and not point_data + Gnum, den = _slot_map(ref_el, entity, J, derivative=False, facet=facet) + tab = _tabulate(fiat_element, 0, ell.points, tabs)[(0,) * sd] + T = tab.reshape(tab.shape[0], sd**ell.rank, npts) + N = numpy.einsum("jpq,qc->pcj", T, ell.weights) + N = N.reshape(sd**ell.rank * sd**ell.rank, -1) + Gprod = reduce(numpy.kron, [Gnum] * ell.rank) + row = _sparse_combination(Gprod.ravel(), N[:, cols], tol, symbolic) + return row, reduce(mul, [den] * ell.rank) + + if ell.order == 0: + # Point values pull back exactly + return None + facet = dim == sd - 1 + Gnum, den = _slot_map(ref_el, entity, J, derivative=True, facet=facet) + direction = ell.pullback(Gnum).direction + tab = _tabulate(fiat_element, ell.order, ell.points, tabs) + N = numpy.stack([tab[alpha] @ ell.weights + for alpha in multiindices(sd, ell.order)]) + row = _sparse_combination(direction, N[:, cols], tol, symbolic) + if not avg and len(ell.points) > 1 and facet: + # The reference weights are measure-intrinsic (integral averages), + # so a plain physical integral carries the physical facet measure. + nu = ref_el.compute_scaled_normal(entity) + Knu = adjugate(J).T @ nu + measure = (ref_el.volume_of_subcomplex(sd - 1, entity) + / (nu @ nu) ** 0.5) * (Knu @ Knu) ** 0.5 + row = row * measure + return row, reduce(mul, [den] * ell.order) diff --git a/finat/piola_mapped.py b/finat/piola_mapped.py index 447237e2..4db254ba 100644 --- a/finat/piola_mapped.py +++ b/finat/piola_mapped.py @@ -98,7 +98,7 @@ def entity_dofs(self): def space_dimension(self): return self._space_dimension - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): sd = self.cell.get_spatial_dimension() bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) diff --git a/finat/walkington.py b/finat/walkington.py index b9d11ae3..dc599bbe 100644 --- a/finat/walkington.py +++ b/finat/walkington.py @@ -24,7 +24,7 @@ def __init__(self, cell, degree=5): reduced_dofs[sd-1][entity] = reduced_dofs[sd-1][entity][:1] self._entity_dofs = reduced_dofs - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): # Jacobian at barycenter sd = self.cell.get_spatial_dimension() top = self.cell.get_topology() diff --git a/finat/wuxu.py b/finat/wuxu.py index 39501566..9b519ad0 100644 --- a/finat/wuxu.py +++ b/finat/wuxu.py @@ -15,7 +15,7 @@ def __init__(self, cell, degree=7): cite("WuXu2019") super().__init__(FIAT.WuXuRobustH3NC(cell)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): return wuxu_transformation(self, coordinate_mapping) @@ -26,7 +26,7 @@ def __init__(self, cell, degree=4): cite("WuXu2019") super().__init__(FIAT.WuXuH3NC(cell)) - def basis_transformation(self, coordinate_mapping): + def _basis_transformation(self, coordinate_mapping): return wuxu_transformation(self, coordinate_mapping) diff --git a/finat/zany.py b/finat/zany.py new file mode 100644 index 00000000..0e8ab539 --- /dev/null +++ b/finat/zany.py @@ -0,0 +1,83 @@ +r"""Family mixins for automatically transformed elements. + +The basis transformation itself is derived generically by +:meth:`finat.physically_mapped.PhysicallyMappedElement.basis_transformation`, +which assembles the generalized Vandermonde matrix :math:`B_{ij} = +n_i(\hat\psi_j\circ F^{-1})` from the FIAT dual basis and inverts it by +sparse block back-substitution over its numerically inferred sparsity +pattern. The classes here only carry the family-specific conventions: + +* :class:`ScalarPhysicallyMappedElement`, for scalar elements with an + affine (identity) pullback -- Morley, Hermite, Argyris, Bell. The + physical facet nodes take their normal component along the unit + physical normal and their tangential components along the mapped + reference tangents; the ``avg`` attribute records whether physical + facet moments are integral averages or plain integrals. + +* :class:`PiolaPhysicallyMappedElement`, for vector- or tensor-valued + elements under the (double) contravariant Piola pullback -- + Mardal-Tai-Winther, Johnson-Mercier, Arnold-Winther, Guzman-Neilan. + The roles of the normal and tangential directions are mirrored: the + scaled facet normal is the cofactor image of the reference one, so + pure normal-component moments are invariant, while scaled tangents + map by the Jacobian. Vertex values of tensor-valued fields are + rescaled by :math:`h^{-2}` following the hand-coded Arnold-Winther + and Hu-Zhang conventions. +""" + +from FIAT.finite_element import FiniteElement + +from finat.physically_mapped import PhysicallyMappedElement + + +class ScalarPhysicallyMappedElement(PhysicallyMappedElement): + """Mixin for scalar elements with an affine (identity) pullback.""" + + def _check_mapping(self, fiat_element: FiniteElement) -> None: + """Reject FIAT elements that are not affinely mapped. + + :arg fiat_element: The FIAT element defined on the reference cell. + :raises NotImplementedError: If the pullback is not affine. + """ + mappings = set(fiat_element.mapping()) + if mappings != {"affine"}: + raise NotImplementedError( + f"{type(self).__name__} expects an affine pullback, not {mappings}.") + + +class PiolaPhysicallyMappedElement(PhysicallyMappedElement): + """Mixin for vector- or tensor-valued elements under the (double) + contravariant Piola pullback.""" + + def _check_mapping(self, fiat_element: FiniteElement) -> None: + """Reject FIAT elements that are not (double) contravariant Piola mapped. + + :arg fiat_element: The FIAT element defined on the reference cell. + :raises NotImplementedError: If the pullback is not supported. + """ + mappings = set(fiat_element.mapping()) + if mappings not in ({"contravariant piola"}, {"double contravariant piola"}): + raise NotImplementedError( + f"{type(self).__name__} expects a (double) contravariant " + f"Piola pullback, not {mappings}.") + + def dof_scale(self, node, dim: int, havg): + r"""Return the conditioning rescaling factor of one physical dof. + + On top of the default derivative-order convention, vertex point + evaluations of tensor-valued (rank-2) fields are redefined with a + factor :math:`h^{-2}`, following the hand-coded Arnold-Winther and + Hu-Zhang transformations; facet moments are already scaled by the + facet measure in FIAT and need no further rescaling. + + :arg node: The FIAT functional of the dof. + :arg dim: Topological dimension of the entity the dof sits on. + :arg havg: GEM scalar for the cell size averaged over the vertices + of the dof's entity. + :returns: The GEM scaling factor, or ``None`` for no rescaling. + """ + if dim == 0 and node.pt_dict: + comps = {comp for pt in node.pt_dict for w, comp in node.pt_dict[pt]} + if len(max(comps)) == 2: + return havg**(-2) + return super().dof_scale(node, dim, havg) diff --git a/gem/AGENTS.md b/gem/AGENTS.md new file mode 100644 index 00000000..7c52b437 --- /dev/null +++ b/gem/AGENTS.md @@ -0,0 +1,29 @@ +### Confusing GEM patterns, clarified + +* **Use operator overloading, not manual node construction.** `gem.Node` overloads + `+ - * / ** @` (via `as_gem`/`componentwise`) to work transparently across GEM nodes, + `gem.Literal`, and plain Python/numpy numbers, with automatic `Zero` folding. Write + `havg**(-m)`, never `gem.Power(havg, gem.Literal(-m))`: the operator form also works when + `havg` is a raw float (as in the test harness's `MyMapping.cell_size()`) whereas a manual + `gem.Power` call assumes GEM operands and will not always coerce correctly. +* **Never call `numpy.linalg.inv`/`solve` on a GEM-valued matrix.** LAPACK cannot see inside + GEM expressions. Symbolic linear solves use `adjugate(A) / determinant(A)` + (`finat/physically_mapped.py`) — the symbolic Cramer's-rule equivalent — e.g. in + `FacetFrame.decompose` and `_piola_facet_rows`. Numeric `numpy.linalg` calls are only valid + when every entry is a plain number: reference-cell-only quantities + (`FacetFrame.reference_coefficients`) or purely numeric direction-basis inversions + (`_scalar_point_rows`, `_piola_point_rows`). Check which regime an array is in before + reaching for a numpy linear-algebra routine. +* **Build sparse GEM-valued arrays with `numpy.full(shape, gem.Zero(), dtype=object)`**, never + `numpy.zeros(shape)` — plain `0` is not interchangeable with `gem.Zero()` when the array will + later be combined with GEM nodes via `+`. +* **"Start from identity, mutate only the rows that need work" is not an optimization, it is + the mathematical content.** `V = identity(ndof)` (`finat/physically_mapped.py`) encodes + "these dofs are push-forward invariant" directly; the invariant-dof detection + (`_invariant_dofs`) simply chooses *not* to touch those rows, rather than writing `1`s + explicitly. Treat the untouched identity rows as the base case of the assembly recursion. +* **Row/column convention, and where the one transpose happens.** Throughout assembly, row + index = reference node, column index = physical node — i.e. the code builds $V$, never $M$ + directly. `ListTensor(V.T)` at the very end is the single place Kirby (2017) Theorem 3.1's + $M = V^T$ gets applied. If something looks transposed, check this convention before + suspecting a sign error. diff --git a/rck_remarks_for_claude.md b/rck_remarks_for_claude.md new file mode 100644 index 00000000..5ec2afd0 --- /dev/null +++ b/rck_remarks_for_claude.md @@ -0,0 +1,53 @@ +We shouldn't assume that det(J) > 0. The UFC ordering we are explicitly utilizing means that we will have both signs. Even in 2d, some triangles will be ordered counterclockwise and others clockwise. This appears in a few places in the document. Note that FIAT's reference normals may be inward or outward pointing (based on UFC convention) and that the existing zany maps in FInAT work without this assumption. + +We need to define "facet frame" explicitly. + + + +Every functional in (2.8) has exactly one direction and order of differentiation. This is probably sufficient for all the nodes we're interested in, but we may find a situation where this needs to be generalized. Don't generalize yet just to do so. + +After (2.9), we say "dividing by J". Change that to multiplying by J^{-T}. + +Be more clear in this paragraph -- push-forward changes the direction and so the actual transformation is geometrically dependent, even if we are doing a lot of preprocessing that can be done numerically. + +Be more clear in the "The design equation" -- we say "conversely", but what are we contrasting with. + +Be more clear about (2.10) -- it's "solvable", but are we going to do that numerically for coefficients of the directions or do we have to do symbolic processing (this will be clear later, but we need to cue the reader as to what's coming) + +Is the direction tensor d always going to be symmetric, and if it's outer(n, n) (say) in 2d, then it has rank 1 instead of 2. Is this reduced rank something that generalizes to other situations (e.g. Wu-Xu) and does it give us anything important to consider? + +Above (2.11) -- "scaled tangent" is FIAT-internal lingo (unit vector times length/measure). Worth defining explicitly. + +Above (2.11) point #3 -- I presume this includes integral moments of derivatives as well? + +Another point in this discussion -- if we're doing things numerically, without dispatching on type, how do we know that we have a normal component.  I think the point is that we don't know and don't need to.  Instead, we're giving suitable bases in reference and physical configurations in which to expand whatever direction we're working with. + +Notation: using big N for the physical normal and little \hat{n} for reference seems odd. Of course, we count things n. Maybe using \nu and \hat{\nu} for physical/reference normals so we don't risk confusion between N and the set of nodes \mathcal{N}? + + +For facet nodes and completion, does this work for derivative nodes that are linear combinations of point derivatives (e.g. integral moments)? +In that case, do we reproduce the fundamental theorem of calculus, writing the completed derivative as the difference between point values on the edges? +This is claimed after (2.21), but a proof that we get FTOC would be nice. + + +Is there a nice analog of some Stokes-like theorem for an integral moment of a normal derivative on a face of a tetrahedron? I think there is something like this in Xu's papers on generalized nonconforming elements. Again, it would be nice to show (a theorem?) that we get that from our numerics in the special case. + +Is it necessary at all to "bin" nodes that have the same set of points and/or weights, or does this sort of fall out numerically? + + + +I also think the recursive reduction to form the V matrix product without E V^c D being all explicit is very nice. However, I think it's worth separating the +presentation into two parts: +1) exactly describing how we can form the three bits E V^c and D separately -- V^c essentially follows from the affine-interpolation equivalent discussion, and D is described by the numerical systems we solve, for example. This could be more explicit. +2) then keep describing the reduction via row recursion as a separate step. It's got to be morally equivalent to some graph algorithm for a sparse triple product, but with blocks based on cell topology somehow. + +(2.23) again may have to deal with the possibility that these matrices have negative determinant? + +Phrase the recursion in pseudocode (algorithmic in latex) + +It would be nice to work out (explicitly with formulas?) the Morley and maybe hexic Argryis to demonstrate the process. +Then, we should put in an example of an element whose transformation is *not* explicitly known in the literature (e.g. one of Xu's generalized nonconforming elements?) + + +Stylistic point: is it worth identifying propositions/theorems about the calcultaions we do to, say, identify the structure/entries of V^c, E, D, and the recursive algorithm for the product? + diff --git a/test/finat/test_claude_zany_piola.py b/test/finat/test_claude_zany_piola.py new file mode 100644 index 00000000..6220d6a6 --- /dev/null +++ b/test/finat/test_claude_zany_piola.py @@ -0,0 +1,466 @@ +r"""Numeric prototype of the unified transformation theory, scalar and Piola. + +The transformation matrix is obtained by duality alone: + +.. math:: B_{ij} = n_i(\hat\psi_j \circ F^{-1}), \qquad V = B^{-1}, + +where :math:`\hat\psi_j` is the reference nodal basis, :math:`F` the cell +map, and :math:`n_i` the physical node. Two ingredients make the rows of +:math:`B` computable without any frame algebra: + +* **Physical nodes by per-slot maps.** The physical node shares the + points and weights of its reference partner; only the directional data + changes, one tensor slot at a time. A derivative slot of a facet node + maps its unit-normal component to the unit physical normal + :math:`K\hat{n}/|K\hat{n}|` (:math:`K = \operatorname{adj}(J)^T` the + cofactor matrix, which maps normals to normals) and its tangential + complement by :math:`J` (mapped tangents); away from facets derivative + slots are Cartesian and keep their reference directions. A + contravariant value slot of a facet moment maps its scaled-normal + component by :math:`K` (the cofactor lemma :math:`K\hat\nu^s = \nu^s` + is exact) and its tangential complement by :math:`J`. Cartesian point + data keeps its weights, interior moments are invariant by convention + (physical test functions are Piola-mapped), and divergence nodes + contract to :math:`\det J` times the identity. + +* **The adjoint acts on the tabulation.** The push-forward of a + reference node contracts each derivative slot of its direction with + :math:`J` (:math:`d = J^{\otimes m}\hat{d}`); dually, instead of + transforming directions, each derivative slot of the *numeric* + reference tabulation is contracted with :math:`J^{-T}` (the physical + derivatives of :math:`\hat\psi_j\circ F^{-1}`) and each value slot + with :math:`J/\det J` (its physical Piola values), once per node + group. A row of :math:`B` is then the plain numeric pairing of the + physical node data with the transformed tabulation; push-forward + invariance (mapped-tangential moments, point values) needs no special + cases -- it falls out as exact Kronecker rows. + +Crucially, :math:`B` is never inverted densely. Its rows obey the support +law :math:`B_{ij} \ne 0` only if dof :math:`j` lives on the closure of dof +:math:`i`'s entity: interior-moment, divergence, and single-point rows are +*exactly* block-diagonal (pointwise identities, independent of the +polynomial space); a vertex-jet row is an exact combination of the +reference jet dofs at the same vertex (duality of the nodal basis); and a +facet row couples to its own facet block plus, possibly, dofs on the +boundary of that facet -- the residual left by the tangential components +is a functional of the trace on the facet, and the trace is determined by +the closure dofs (the same property that makes the element conforming; in +the scalar theory this role is played by the fundamental-theorem-of- +calculus exactness identities). Hence :math:`B` is block lower triangular +in the entity partial order and :func:`composition_transformation` +computes :math:`V = B^{-1}` by block back-substitution -- one small solve +per entity, fill-in confined to the closure. Both the support law and the +agreement with ``basis_transformation`` are asserted to machine precision +for the scalar and Piola element zoos, in 2D and 3D, on cells of both +orientations. This module is the independent verification of the +automated transformations: it shares no code path with +``basis_transformation``. See ``zany_claude.md`` (Stages 4-5) for the +derivation. +""" + +from math import factorial, prod + +import numpy as np +import pytest + +import finat +from FIAT.polynomial_set import mis +from FIAT.reference_element import make_affine_mapping, ufc_simplex +from finat.functional import PhysicallyMappedFunctional +from gem.interpreter import evaluate + +from .conftest import MyMapping + + +def contract_slot(T, A, axis): + """Contract one tensor slot with a matrix. + + :arg T: The tensor. + :arg A: The matrix. + :arg axis: The slot of ``T`` to contract. + :returns: The tensor with ``T'[..., i, ...] = sum_k A[i, k] T[..., k, ...]``. + """ + return np.moveaxis(np.tensordot(T, A, axes=(axis, 1)), -1, axis) + + +def tensorize_direction(direction, sd, order): + """Distribute multi-index direction coefficients over a symmetric tensor. + + :arg direction: Direction coefficients in the multi-index basis. + :arg sd: The spatial dimension. + :arg order: The derivative order. + :returns: The symmetric direction tensor, normalized so that the full + contraction with an index-wise derivative tabulation reproduces + the multi-index pairing. + """ + alphas = sorted(mis(sd, order), reverse=True) + lookup = {alpha: k for k, alpha in enumerate(alphas)} + D = np.zeros((sd,) * order) + for index in np.ndindex(D.shape): + alpha = [0] * sd + for k in index: + alpha[k] += 1 + scale = prod(map(factorial, alpha)) / factorial(order) + D[index] = direction[lookup[tuple(alpha)]] * scale + return D + + +def derivative_tabulation(fiat_element, order, points, weights): + """Weighted symmetric-tensor derivative tabulation of the nodal basis. + + :arg fiat_element: The FIAT element. + :arg order: The derivative order. + :arg points: The quadrature points. + :arg weights: The quadrature weights. + :returns: The tensor ``T[j, i1, ..., im] = sum_q w_q + (d^m psi_j / dx_i1 ... dx_im)(x_q)``. + """ + sd = fiat_element.get_reference_element().get_spatial_dimension() + tab = fiat_element.tabulate(order, points) + nbf = tab[(0,) * sd].shape[0] + T = np.zeros((nbf,) + (sd,) * order) + for index in np.ndindex((sd,) * order): + alpha = [0] * sd + for k in index: + alpha[k] += 1 + T[(slice(None), *index)] = tab[tuple(alpha)] @ weights + return T + + +def facet_slot_map(ref_el, entity, J, derivative): + r"""Per-slot physical direction map of a facet node. + + A derivative slot maps its component along the FIAT normal + :math:`\hat{n}` to the FIAT physical normal -- the norm-preserving + rescaling of the cofactor image :math:`K\hat{n}` (:math:`K = + \operatorname{adj}(J)^T` maps normals to normals, and the norm of + the FIAT normal depends only on the reference cell) -- and its + tangential complement by :math:`J` (mapped tangents). + + A contravariant value slot maps its component along the scaled + normal :math:`\hat\nu^s` by :math:`K` (the cofactor lemma + :math:`K\hat\nu^s = \nu^s` is exact), and its tangential complement + to the cofactor image projected onto the physical facet, scaled by + :math:`|K\hat\nu^s|^2/(\det J\, |\hat\nu^s|^2)` -- the reciprocal + of the scalar tangential push-forward coefficient, which in 2D + reduces to the mapped tangent :math:`J\hat{t}`. + + :arg ref_el: The reference cell. + :arg entity: The facet number. + :arg J: The (numeric) cell Jacobian. + :arg derivative: If True, use the derivative-slot convention; + if False, the contravariant value-slot convention. + :returns: The map as an ``(sd, sd)`` array acting on one slot. + """ + sd = ref_el.get_spatial_dimension() + detJ = np.linalg.det(J) + K = detJ * np.linalg.inv(J).T + if derivative: + n = ref_el.compute_normal(entity) + Kn = K @ n + Kn = Kn * (np.linalg.norm(n) / np.linalg.norm(Kn)) + P = np.outer(n, n) / (n @ n) + return np.outer(Kn, n) / (n @ n) + J @ (np.eye(sd) - P) + n = ref_el.compute_scaled_normal(entity) + Kn = K @ n + P = np.outer(n, n) / (n @ n) + Q = np.eye(sd) - np.outer(Kn, Kn) / (Kn @ Kn) + s = (Kn @ Kn) / (detJ * (n @ n)) + return np.outer(Kn, n) / (n @ n) + s * (Q @ K @ (np.eye(sd) - P)) + + +def scalar_derivative_row(fiat_element, ell, Phi, J): + r"""B row of a scalar derivative node. + + The physical direction is the per-slot map :math:`\Phi^{\otimes m}` + of the reference direction; the derivative slots of the tabulation + carry the adjoint of the push-forward, :math:`J^{-T}` per slot. + + :arg fiat_element: The FIAT element. + :arg ell: The parsed reference node. + :arg Phi: The per-slot physical direction map. + :arg J: The (numeric) cell Jacobian. + :returns: The vector of values of the physical node on the + pushed-forward nodal basis. + """ + sd = J.shape[0] + T = derivative_tabulation(fiat_element, ell.order, ell.points, ell.weights) + A = np.linalg.inv(J).T + D = tensorize_direction(ell.direction, sd, ell.order) + for slot in range(ell.order): + T = contract_slot(T, A, 1 + slot) + D = contract_slot(D, Phi, slot) + return np.tensordot(T, D, axes=(tuple(range(1, ell.order + 1)), + tuple(range(ell.order)))) + + +def piola_value_row(fiat_element, ell, Phi, J): + r"""B row of a contravariant value node. + + The value slots of the tabulation carry the physical Piola values, + :math:`J/\det J` per slot; the weights are mapped per slot by + ``Phi`` (facet moments) or kept Cartesian (point data). + + :arg fiat_element: The FIAT element. + :arg ell: The parsed reference node. + :arg Phi: The per-slot physical weight map, or None for Cartesian + point data. + :arg J: The (numeric) cell Jacobian. + :returns: The vector of values of the physical node on the + pushed-forward nodal basis. + """ + sd = J.shape[0] + r = ell.rank + npts = len(ell.points) + Theta = J / np.linalg.det(J) + T = fiat_element.tabulate(0, ell.points)[(0,) * sd] + T = T.reshape(T.shape[0], *(sd,) * r, npts) + W = ell.weights.reshape(npts, *(sd,) * r) + for slot in range(r): + T = contract_slot(T, Theta, 1 + slot) + if Phi is not None: + W = contract_slot(W, Phi, 1 + slot) + return np.tensordot(T, W, axes=((*range(1, r + 1), r + 1), + (*range(1, r + 1), 0))) + + +def evaluate_divergence(ell, fiat_element): + """Apply a divergence functional to the nodal basis of a FIAT element. + + :arg ell: A divergence :class:`PhysicallyMappedFunctional`. + :arg fiat_element: The FIAT element providing the nodal basis. + :returns: The vector of values of the functional on the nodal basis. + """ + sd = fiat_element.get_reference_element().get_spatial_dimension() + tab = fiat_element.tabulate(1, ell.points) + alphas = [tuple(int(k == c) for k in range(sd)) for c in range(sd)] + div = sum(tab[alphas[c]].reshape(tab[alphas[c]].shape[0], sd, -1)[:, c, :] + for c in range(sd)) + return div @ ell.weights + + +def physical_node_row(fiat_element, ell, dim, entity, J, avg): + """B row of one parseable node, or None if push-forward invariant. + + :arg fiat_element: The FIAT element. + :arg ell: The parsed reference node. + :arg dim: The dimension of the entity the node sits on. + :arg entity: The entity number. + :arg J: The (numeric) cell Jacobian. + :arg avg: If False, physical scalar facet moments are plain + integrals rather than the measure-intrinsic integral averages of + the reference weights, and the row is rescaled by the physical + facet measure. + :returns: The row, or None for a row of the identity. + """ + ref_el = fiat_element.get_reference_element() + sd = ref_el.get_spatial_dimension() + if ell.divergence: + return evaluate_divergence(ell, fiat_element) / np.linalg.det(J) + if ell.rank: + point_data = len(ell.points) == 1 and (ell.rank == 1 or dim != sd - 1) + if dim == sd and not point_data: + return None + Phi = (facet_slot_map(ref_el, entity, J, derivative=False) + if dim == sd - 1 and not point_data else None) + return piola_value_row(fiat_element, ell, Phi, J) + if ell.order == 0: + return None + Phi = (facet_slot_map(ref_el, entity, J, derivative=True) + if dim == sd - 1 else np.eye(sd)) + row = scalar_derivative_row(fiat_element, ell, Phi, J) + if not avg and len(ell.points) > 1 and dim == sd - 1: + # The reference weights are measure-intrinsic (integral averages), + # so a plain physical integral carries the physical facet measure. + n = ref_el.compute_scaled_normal(entity) + K = np.linalg.det(J) * np.linalg.inv(J).T + measure = (ref_el.volume_of_subcomplex(sd - 1, entity) + * np.linalg.norm(K @ n) / np.linalg.norm(n)) + row = row * measure + return row + + +def closure_dofs(fiat_element, dim, entity): + """Indices of the dofs on the strict closure of a cell entity. + + :arg fiat_element: The FIAT element. + :arg dim: The dimension of the entity. + :arg entity: The entity number. + :returns: The indices of the dofs supported on subentities of strictly + lower dimension contained in the closure of the entity. + """ + top = fiat_element.get_reference_element().get_topology() + entity_ids = fiat_element.entity_dofs() + verts = set(top[dim][entity]) + return [i for d in sorted(entity_ids) if d < dim + for e in entity_ids[d] if set(top[d][e]) <= verts + for i in entity_ids[d][e]] + + +def composition_transformation(fiat_element, J, ndof=None, avg=True, tol=2e-10): + """Compute V for a scalar or Piola element by blockwise duality. + + The matrix :math:`B` of the module docstring is block lower triangular + in the entity partial order, so :math:`V = B^{-1}` is computed by block + back-substitution over the entities in order of increasing dimension: + + .. math:: V_e = B_{ee}^{-1} (I_e - B_{ec} V_c), + + where :math:`e` collects the dofs of one entity and :math:`c` the + (already processed) dofs on its strict closure. The support law that + justifies the triangular structure is asserted, not assumed: every row + of :math:`B` must vanish outside its own entity block and closure. + Fill-in is confined to the closure, so the sparsity of :math:`V` + matches that of the row recursion in the scalar theory. + + :arg fiat_element: The FIAT element on the reference cell. + :arg J: The (numeric) cell Jacobian. + :arg ndof: The number of exposed physical dofs; unparseable constraint + functionals beyond it keep identity rows (their columns are + truncated by the caller). Defaults to all dofs. + :arg avg: Whether physical scalar facet moments are integral averages. + :arg tol: Relative tolerance for the support-law assertion. + :returns: The transformation V as a numpy array, with the same row and + column ordering as ``basis_transformation`` before truncation of + the trailing constraint columns. + """ + nodes = fiat_element.dual_basis() + entity_ids = fiat_element.entity_dofs() + nbf = len(nodes) + if ndof is None: + ndof = nbf + eye = np.eye(nbf) + V = np.zeros((nbf, nbf)) + done = set() + for dim in sorted(entity_ids): + for entity in entity_ids[dim]: + # FIAT may list a dof on more than one entity (e.g. the edge + # moments of Arnold-Winther reappear in its interior list); the + # lowest-dimensional entity owns the dof. + block = [i for i in entity_ids[dim][entity] if i not in done] + if not block: + continue + done.update(block) + rows = [] + for i in block: + try: + ell = PhysicallyMappedFunctional.from_fiat(nodes[i]) + except NotImplementedError: + if i < ndof: + raise + rows.append(eye[i]) + continue + row = physical_node_row(fiat_element, ell, dim, entity, J, avg) + rows.append(eye[i] if row is None else row) + Bblock = np.asarray(rows) + prior = closure_dofs(fiat_element, dim, entity) + outside = np.setdiff1d(np.arange(nbf), block + prior) + if outside.size: + assert (np.abs(Bblock[:, outside]).max() + <= tol * np.abs(Bblock).max()), \ + f"support law violated on entity ({dim}, {entity})" + V[block] = np.linalg.solve( + Bblock[:, block], eye[block] - Bblock[:, prior] @ V[prior]) + return V + + +# Scalar elements, including hand-coded macroelement transformations +# (HCT, Powell-Sabin, Alfeld/Bramble-Zlamal C2) never before reproduced by +# the automated theory. BrambleZlamalC2 needs a looser tolerance: FIAT's +# own dual basis and tabulation only agree to ~4e-10 on its ill-conditioned +# order-4 vertex jets. Walkington is excluded: its extended tabulation has +# full rank, so the constraint columns of the transformation are not a +# convention, and the hand-coded transformation disagrees there both with +# this prototype and with a direct physical-cell fit (see zany_claude.md). +scalar_zoo = { + 2: [(finat.Morley, ()), + (finat.Hermite, ()), + (finat.Bell, ()), + (finat.Argyris, ()), + (finat.Argyris, (5, "point")), + (finat.Argyris, (6,)), + (finat.WuXuH3NC, ()), + (finat.WuXuRobustH3NC, ()), + (finat.HsiehCloughTocher, ()), + (finat.HsiehCloughTocher, (4,)), + (finat.ReducedHsiehCloughTocher, ()), + (finat.QuadraticPowellSabin6, ()), + (finat.QuadraticPowellSabin12, ()), + (finat.AlfeldC2, ()), + (finat.BrambleZlamalC2, ()), + ], + 3: [(finat.Morley, ()), + (finat.Hermite, ())], +} + +piola_zoo = { + 2: [(finat.MardalTaiWinther, ()), + (finat.JohnsonMercier, ()), + (finat.ArnoldWintherNC, ()), + (finat.ArnoldWinther, ()), + (finat.AlfeldSorokina, ()), + (finat.BernardiRaugel, ()), + (finat.BernardiRaugelBubble, ()), + (finat.GuzmanNeilanFirstKindH1, ()), + (finat.GuzmanNeilanSecondKindH1, ()), + (finat.GuzmanNeilanBubble, ()), + (finat.GuzmanNeilanH1div, ()), + (finat.ReducedArnoldQin, ()), + (finat.ChristiansenHu, ()), + (finat.HuZhang, (3, "integral")), + (finat.HuZhang, (4, "integral")), + (finat.HuZhang, (3, "point")), + (finat.HuZhang, (4, "point"))], + 3: [(finat.MardalTaiWinther, ()), + (finat.MardalTaiWinther, (2,)), + (finat.JohnsonMercier, ()), + (finat.AlfeldSorokina, ()), + (finat.BernardiRaugel, ()), + (finat.BernardiRaugelBubble, ()), + (finat.GuzmanNeilanFirstKindH1, ()), + (finat.GuzmanNeilanFirstKindH1, (2,)), + (finat.GuzmanNeilanSecondKindH1, ()), + (finat.GuzmanNeilanBubble, ()), + (finat.GuzmanNeilanH1div, ()), + (finat.ChristiansenHu, ())], +} + +orientations = { + (2, "positive"): ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84)), + (2, "negative"): ((0.0, 0.1), (0.15, 1.84), (1.17, -0.09)), + (3, "positive"): ((0, 0, 0), (1., 0.1, -0.37), + (0.01, 0.987, -.23), (-0.1, -0.2, 1.38)), + (3, "negative"): ((0, 0, 0), (0.01, 0.987, -.23), + (1., 0.1, -0.37), (-0.1, -0.2, 1.38)), +} + + +def check_composition(dimension, element, args, orientation): + """Compare the blockwise duality V against basis_transformation.""" + ref_cell = ufc_simplex(dimension) + phys_cell = ufc_simplex(dimension) + phys_cell.vertices = orientations[(dimension, orientation)] + J, b = make_affine_mapping(ref_cell.vertices, phys_cell.vertices) + + finat_element = element(ref_cell, *args) + mapping = MyMapping(ref_cell, phys_cell) + M = evaluate([finat_element.basis_transformation(mapping)])[0].arr + + ndof = finat_element.space_dimension() + avg = getattr(finat_element, "avg", True) + V = composition_transformation(finat_element._element, J, ndof=ndof, avg=avg) + assert np.allclose(V[:, :ndof], M.T, atol=2e-10) + + +@pytest.mark.parametrize("orientation", ["positive", "negative"]) +@pytest.mark.parametrize("dimension, element, args", [ + (dim, *case) for dim in scalar_zoo for case in scalar_zoo[dim]]) +def test_scalar_composition(dimension, element, args, orientation): + check_composition(dimension, element, args, orientation) + + +@pytest.mark.parametrize("orientation", ["positive", "negative"]) +@pytest.mark.parametrize("dimension, element, args", [ + (dim, *case) for dim in piola_zoo for case in piola_zoo[dim]]) +def test_piola_composition(dimension, element, args, orientation): + check_composition(dimension, element, args, orientation) diff --git a/test/finat/test_zany_automation.py b/test/finat/test_zany_automation.py new file mode 100644 index 00000000..e68f74ca --- /dev/null +++ b/test/finat/test_zany_automation.py @@ -0,0 +1,79 @@ +import FIAT +import finat +import numpy as np +import pytest + +from gem.interpreter import evaluate +from finat.functional import PhysicallyMappedFunctional +from finat.fiat_elements import FiatElement +from finat.zany import PiolaPhysicallyMappedElement, ScalarPhysicallyMappedElement + + +@pytest.mark.parametrize("dimension", [2, 3]) +def test_functional_from_fiat(dimension): + """The symbolic PhysicallyMappedFunctional recovers order, weights and + direction numerically from the FIAT functional dictionaries.""" + cell = FIAT.ufc_simplex(dimension) + element = FIAT.Morley(cell) + entity_ids = element.entity_dofs() + nodes = element.dual_basis() + + for i in entity_ids[dimension - 2][0]: + ell = PhysicallyMappedFunctional.from_fiat(nodes[i]) + assert ell.order == 0 + assert ell.direction is None + + for entity in entity_ids[dimension - 1]: + for i in entity_ids[dimension - 1][entity]: + ell = PhysicallyMappedFunctional.from_fiat(nodes[i]) + assert ell.order == 1 + normal = cell.compute_normal(entity) + cosine = ell.direction @ normal + assert np.isclose(abs(cosine), np.linalg.norm(normal)) + + +auto_elements = [finat.Morley, finat.Hermite] + + +@pytest.mark.parametrize("element", auto_elements) +@pytest.mark.parametrize("dimension", [2, 3]) +def test_conditioning_scaling(ref_to_phys, scaled_ref_to_phys, element, dimension): + """Derivative dofs are rescaled by cell size to the power of the + derivative order.""" + scaled = scaled_ref_to_phys[dimension][-1] + # the same geometric mapping with unit cell size + unit = type(ref_to_phys[dimension])(scaled.ref_cell, scaled.phys_cell) + finat_element = element(scaled.ref_cell) + + Ms = evaluate([finat_element.basis_transformation(scaled)])[0].arr + Mu = evaluate([finat_element.basis_transformation(unit)])[0].arr + + h = scaled.cell_size()[0] + assert not np.isclose(h, 1) + orders = [node.max_deriv_order for node in finat_element._element.dual_basis()] + expected = Mu * np.asarray([h**-order for order in orders])[:, None] + assert np.allclose(Ms, expected) + + +def test_unsupported_mapping_scalar(ref_to_phys): + """A ScalarPhysicallyMappedElement rejects a non-affine FIAT element.""" + mapping = ref_to_phys[2] + + class BogusScalar(ScalarPhysicallyMappedElement, FiatElement): + pass + + element = BogusScalar(FIAT.Nedelec(mapping.ref_cell, 1)) + with pytest.raises(NotImplementedError): + element.basis_transformation(mapping) + + +def test_unsupported_mapping_piola(ref_to_phys): + """A PiolaPhysicallyMappedElement rejects a covariant FIAT element.""" + mapping = ref_to_phys[2] + + class BogusPiola(PiolaPhysicallyMappedElement, FiatElement): + pass + + element = BogusPiola(FIAT.Nedelec(mapping.ref_cell, 1)) + with pytest.raises(NotImplementedError): + element.basis_transformation(mapping) diff --git a/zany_claude.md b/zany_claude.md new file mode 100644 index 00000000..b79eb3e0 --- /dev/null +++ b/zany_claude.md @@ -0,0 +1,847 @@ +## Transformation Theory (Kirby 2017; Brubeck & Kirby 2025) + +The mathematical framework for mapping non-affinely-equivalent elements. Notation: the +geometric map goes **from physical to reference**, $F: K \to \hat{K}$; pullback +$F^*(\hat{f}) = \hat{f}\circ F$; push-forward of a node $F_*(n) = n \circ F^*$. + +* Goal: the matrix $M$ with $\Psi = M\, F^*(\hat\Psi)$, expressing physical nodal basis + functions as combinations of pulled-back reference basis functions. +* Duality is the whole game: with $B_{ij} = F_*(n_i)(\hat\psi_j) = n_i(F^*(\hat\psi_j))$, + the Kronecker property gives $M = B^{-T}$ and $V = B^{-1}$, where $V$ relates nodes: + $\hat{N} = V\, F_*(N)$ (restricted to $P$). Hence **$M = V^T$**, and $V$ is what one + actually constructs, because push-forwards of nodes are computable by the chain rule. +* Affine equivalence (Lagrange): $F_*(N) = \hat{N}$, so $M = I$. +* Affine-interpolation equivalence (Hermite): spans of nodes at each point are preserved, + so $V$ is block diagonal, one small block per point-node group (e.g. $J^{-T}$ blocks for + vertex gradients — in the *paper's* $J$; see the FInAT convention below). +* Neither (Morley, Argyris, HCT): edge normal derivatives alone do not push forward into + the span of reference nodes. Fix by a **compatible nodal completion** $N^c \supset N$, + $\hat{N}^c \supset \hat{N}$ with $\mathrm{span}(F_*(N^c)) = \mathrm{span}(\hat{N}^c)$ + (add the tangential-derivative partners so each edge carries a full gradient). Then + $$V = E\, V^c\, D,$$ + with the three factors having distinct mathematical roles: + * $D \in \mathbb{R}^{\mu\times\nu}$: expresses the completed physical nodes in terms of + the actual physical nodes, *restricted to* $P$ ($\pi N^c = D\, \pi N$). Rows for nodes + already in $N$ are Boolean. Rows for completion nodes come from a univariate exactness + argument on the edge: a tangential-derivative quantity of a polynomial of known edge + degree is an exact linear combination of the endpoint values/derivatives (FTC for HCT's + $\mu^t_e$; the quintic endpoint rule for Argyris/Bell midpoint tangential derivatives). + This is the only factor that uses the polynomial degree of the space. + * $V^c \in \mathbb{R}^{\mu\times\mu}$: block diagonal, pure chain rule, + $\hat{N}^c = V^c F_*(N^c)$. Vertex value $\to 1$; vertex gradient $\to$ Jacobian block; + vertex Hessian $\to$ symmetric-square block $\Theta$; edge normal/tangential pair $\to$ + the $2\times2$ block $B_i = \hat{G}_i J^{-T} G_i^T$ (paper's $J$), where + $G = [\mathbf{n}\; \mathbf{t}]^T$. Only the first row of $B_i$ is needed + ($B_{nn}$, $B_{nt}$) since tangential rows are discarded by $E$. + * $E \in \mathbb{R}^{\nu\times\mu}$: Boolean extraction of $\hat{N}$ from $\hat{N}^c$. +* Constrained spaces, $F^*(\hat{P}) \neq P$ (Bell, reduced HCT): the space is + $P = \bigcap_i \mathrm{null}(\lambda_i)$ inside a larger $\tilde{P}$ that *is* preserved + (e.g. reduced HCT = HCT functions whose edge normal derivative is linear; + $\lambda_i$ = moment of the normal derivative against the quadratic Legendre polynomial + on edge $i$). Build the **extended element** $(K, \tilde{P}, [N; L])$ with the + constraints $\lambda_i$ appended as extra nodes; it is a valid finite element, and its + first $\nu$ nodal basis functions are exactly the nodal basis of the constrained + element. Transform the extended element with the $E V^c D$ machinery (completing each + $\lambda_i$ with its tangential partner $\lambda_i'$, again eliminated through an edge + exactness rule), then discard the constraint rows. +* Brubeck & Kirby 2025 refinements: + * Define *reference* edge nodes as integral **averages** ($1/|\hat{e}|$ scaling) while + physical nodes are plain moments: the edge-length Jacobian of the line integral is then + absorbed and no reference-edge-length bookkeeping or orientation logic is needed. + * High-order HCT/Argyris: edge normal-derivative moments are taken against Jacobi + polynomials $P_i^{(1,1)}$ (weights $(2,2)$ for Argyris, matching the vertex jet order), + so the edge-based functions are hierarchical. The tangential completion + $\mu^t_{e,i}$ integrates by parts against the *trace* moments (defined against + $\frac{d}{ds}P_i^{(1,1)}$): $\mu^t_{e,i}(f) = -\mu_{e,i}(f) + P_i^{(1,1)}(1)\, + \delta_{v_b}(f) - P_i^{(1,1)}(-1)\,\delta_{v_a}(f)$. This couples normal moments to + trace moments, not just vertex dofs. + +### Active work: automating the transformation (see `plan_zany_auto.md`) + +Goal: replace the hand-coded `basis_transformation` methods with a helper that derives +$V = E V^c D$ directly from a FIAT element's dual basis. The implementation must mirror +the theory factor by factor, not merely reproduce matrix entries: + +**Status (2026-07-13).** The symbolic dof lives in `finat/functional.py` as +`finat.PhysicallyMappedFunctional`. `finat/physically_mapped.py` stays fully generic: +`PhysicallyMappedElement` is unchanged from before this project, still an abstract +mixin with no knowledge of the zany theory, used as-is by hand-coded elements (AW, +HCT, PowellSabin, Walkington, ...). All the automation lives in `finat/zany.py`, as a +template method on `ZanyPhysicallyMappedElement(PhysicallyMappedElement)`: the +entity-by-entity assembly loop is implemented once, in its concrete +`basis_transformation`, calling four hooks (`_check_mapping`, `_invariant_dofs`, +`_facet_dof_rows`, `_point_dof_rows`) that carry ALL mapping-specific knowledge — the +loop itself contains no `if piola` anywhere. Two mixins implement those hooks, +`ScalarPhysicallyMappedElement` (affine pullback: Morley, Hermite, Argyris, Bell) and +`PiolaPhysicallyMappedElement` ((double) contravariant Piola: MTW, Johnson-Mercier, +Guzman-Neilan), plus the pure math functions they call (`FacetFrame`, +`_scalar_facet_rows`, `_scalar_point_rows`, `_piola_facet_rows`, `_piola_point_rows`) — +these take plain arrays/GEM expressions, no `self`, so the mathematics stays readable +independent of the class plumbing. Concrete elements (`finat.Morley`, etc.) are now +just a citation plus a FIAT constructor call: mixing in the right base class is enough, +`basis_transformation` is inherited (MRO example: `Morley -> ScalarPhysicallyMappedElement +-> ZanyPhysicallyMappedElement -> PhysicallyMappedElement -> ...`). `ndof` truncation is +no longer a parameter; the loop always slices by `self.space_dimension()`, which +constrained elements (Bell, GN) already override. Tests: +`test/finat/test_zany_automation.py`; `check_zany_mapping` lives in the finat conftest +and is provided to test modules as a pytest fixture (pytest runs with +`--import-mode=importlib`, so test modules cannot import from each other or from +conftest). + +**Framework design.** A dof is a symbolic `finat.PhysicallyMappedFunctional`: +$\ell(f) = \sum_q w_q \langle D, \nabla^m f(x_q)\rangle$ with numeric points/weights +and a direction tensor $D$ that is numeric on the reference cell and GEM otherwise. +There is *no dispatch over FIAT functional types*: `PhysicallyMappedFunctional.from_fiat` +reads only `pt_dict`/`deriv_dict` and recovers the order and common direction numerically +(rank-one SVD of the derivative weights). Operations: covariant `pullback(J)` +(contract direction slots with $J$), `with_direction`, and numeric `evaluate` against +a nodal basis (the generalized Vandermonde realizing the $D$ factor of $V = E V^c D$). +The row of $V$ for a reference node $\hat\ell$ with direction $\hat d = a\hat n + +\sum_k \beta_k \hat t_k$ (numeric split) is assembled generically: + +* $a = 0$ (value dofs, or tangential directions): push-forward invariant, identity row. +* otherwise solve $J\hat d = x_0 C + \sum_k x_k J\hat t_k$ symbolically + (`FacetFrame.decompose`, adjugate/determinant of the polynomial frame matrix), where + $C$ = generalized cross product of the mapped tangents. The physical node has + direction $aN + \sum\beta_k J\hat t_k$ with $N = \kappa C/\|C\|$, so its coefficient + is $c = x_0\|C\|/(\kappa a)$ and the tangential remainders are $r_k = x_k - c\beta_k$ + (careful: *not* $a x_k$ — the SVD direction may be $-\hat n$, and all formulas must be + invariant under that sign flip; a sign bug here flips exactly the edges where the SVD + chose the opposite orientation). +* the remainders multiply completion functionals along *mapped reference tangents*, + which coincide with reference functionals; `PhysicallyMappedFunctional.evaluate` gives their numeric + expansion in the element's own nodes, and the row combination recurses through the + 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 +gradients) have no geometric frame: FIAT keeps Cartesian directions on the physical +cell, so the group of derivative nodes on the entity acts as its own completion — this +is precisely affine-interpolation equivalence. The pulled-back direction $J\hat d_i$ +is expanded in the group's own (numeric) direction basis, with weight-ratio factors +making the expansion invariant under the SVD scale/sign ambiguity of each node's +recovered $(w, D)$ factorization. The group must span the derivative jet, share its +points, and have pairwise-parallel weights; otherwise `NotImplementedError`. + +Key facts the framework rests on: + +* **FIAT normals are "UFC consistent":** computed from the tangents by the same formula + on reference and physical cells (`UFCTriangle.compute_normal` rotates the scaled edge + tangent; `UFCTetrahedron.compute_normal` is $-2\times$ the unit cross product of the + scaled face tangents), with cell-independent magnitude. Hence $N = \kappa C/\|C\|$ + with $\kappa$ recoverable from reference data, and no orientation logic is needed + (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 + (`FacetQuadratureRule(..., avg=True)`), so physical nodes share the reference weights + and no `physical_edge_lengths` appear. The framework *assumes* measure-intrinsic + moments (the Brubeck & Kirby 2025 reference-node convention) — documented in + `finat/functional.py`. +* **The conditioning $h$-scaling generalizes via `max_deriv_order`.** Column $j$ is + scaled by $h_E^{-m}$ where $m$ is the derivative order of node $j$ and $h_E$ averages + `cell_size()` over the vertices of its entity; reproduces every per-element + hand-written scaling. `cell_size()` returns raw numpy values in the test mappings but + GEM in Firedrake, so use operator arithmetic (`havg**(-m)`), not GEM constructors + (GEM overloads `+ - * / ** @` with `Zero`/constant folding; keep numpy object arrays + on the left when scaling by a GEM scalar). + +Extensions beyond first order and Morley/Hermite: + +* `PhysicallyMappedFunctional` directions live in derivative multi-index space (`multiindices`, axis + order for $m=1$); `pullback` distributes them over a symmetric tensor (dividing by + multiplicities), contracts every slot with $J$ (`numpy.tensordot` on object arrays), + and collapses back. Point-jet groups are split per order; each order solves in its + own multi-index direction basis (Argyris/Bell vertex jets: gradient + Hessian). +* Facet completions of Argyris edge moments couple to vertex jets and same-edge trace + moments; the existing row recursion handles both with no new code (trace moments are + order-0 and thus invariant; FIAT builds all these moments with + `FacetQuadratureRule(avg=True)`, i.e. measure-intrinsic, as the framework assumes). +* `ScalarPhysicallyMappedElement.avg = False` (an instance attribute Argyris sets from + its constructor kwarg) reproduces the legacy FInAT convention where physical facet + moments are plain integrals: their columns are divided by the physical facet measure + $\|C\| |\hat e| / \|\hat C\|$ (`FacetFrame.measure`). Single-point facet dofs + (Argyris "point" variant) are unaffected. +* Bell is the extended-element pattern: FIAT.Bell is the 21-node quintic element with + the constraint functionals as extra edge nodes; overriding `space_dimension()` to 18 + drops the constraint *columns* of $V$ at the end of the template method (their rows + still contribute the $D$-matrix entries through the completion recursion), and the + FInAT element overrides `entity_dofs`. +* Known convention change: the generic $h^{-m}$ conditioning scaling now also applies + to integral-variant Argyris edge moments, which the hand-written code left unscaled + (Morley scaled them; the legacy convention was inconsistent). Invisible when + `cell_size == 1`; flag in PR review. + +**Piola-mapped elements** (Aznaran, Kirby & Farrell 2022). `PhysicallyMappedFunctional` carries a +value rank: component weight profiles (nq x sd^rank) parsed from `pt_dict` component +tuples. Under contravariant Piola the roles of the scalar case are mirrored: the +*scaled* facet normal is the cofactor image $K\hat n_s$, $K = \mathrm{adj}(J)^T$ +(exactly the physical `compute_scaled_normal`, cross product of mapped tangents), so +pure normal moments are invariant, while scaled tangents map by $J$. `_piola_facet_rows` +works with per-point *frame-coordinate profiles* (handles 3D MTW's point-varying +RT-mapped tangential directions): the pulled-back profile is contracted per value slot +with the mixing matrix $Y$; tangential profiles are matched within the facet group by +solving the (small, square) Gram system $B B^T c = B\cdot(\text{target})$, where $B$ +stacks the group's own reference tangential profiles and has full row rank by +unisolvence, so $B B^T$ is invertible and a rank deficiency (a genuine bug) surfaces as +a hard numerical error rather than passing silently; the residual normal profile is +eliminated by per-point normal moments through the Vandermonde recursion (this is where +e.g. tangential-to-normal couplings emerge). + +**Sparsity gotcha: fold the quadrature-point sum numerically, never symbolically.** The +naive way to do that last elimination is to build the per-point residual (a GEM +expression, since it involves the symbolic mixing matrix $Y$) and dot it against the +tabulated-basis row `L[m, :]` (numeric) with a Python loop over points, accumulating +into `V[i] += V[m] * (residual[q] * L[m, q])`. This reproduces the right *numbers* but +not the right *sparsity*: whenever the true coupling to some dof `m` is zero, it is zero +only as an analytic identity in $J$ (a cancellation across quadrature points or across +frame directions), not as a literal `gem.Zero()` node — GEM has no polynomial-identity +simplifier, so `isinstance(x, gem.Zero)` (or eyeballing `x != 0`) cannot detect this, and +the resulting matrix keeps a dense cloud of numerically-tiny-but-symbolically-nonzero +entries (compare against the hand-coded matrix's exact sparsity, e.g. MTW's +`normal_tangential_transform`, to see the gap). The fix: contract the quadrature-point +axis while everything is still numeric — `Lmap[i] = L @ coords[i]` (both plain arrays) +— and only *afterwards* multiply by the (few, small) symbolic frame-mixing scalars +($Y[0, r]$, or $-c_j$ for each group member's own normal profile), one reference-frame +multi-index/group-member at a time rather than summed together first. Sparsity is then +decided by thresholding the numeric `Lmap` columns (`abs(...) < tol`), exactly the same +idiom already used for `Binv`/`B` above; a symbolic sum of several *different* GEM +scalars must never be built first and then inspected for zero-ness, since that sum's +GEM shape does not reveal whether it is analytically zero. Verified by comparing bit-for- +bit against the deleted hand-coded `MardalTaiWinther.basis_transformation` sparsity +pattern for order 1 (2D/3D) and order 2 (3D): zero extra and zero missing entries. + +Key subtlety (in any dimension > 2): FIAT builds tangential value +components on the **reciprocal basis** (`cross(n, t_k)`), which transforms in-plane +contravariantly: absorb $S^{-1} = (\det\hat G_t/\det G_t)\hat G_t^{-1} G_t$ +(tangent Gram change) into $Y$'s tangential rows; in 2D $S = 1$, which hides the effect. +Interior value moments are Piola-invariant by construction (scaled-normal components for +JM; covariant Nedelec test fields for MTW order 2 cancel the contravariant trial +exactly). Double contravariant (tensor) elements use the same code: one contraction per +value slot. MTW (2D/3D) and Johnson-Mercier (2D/3D) are reimplemented this way, matching +the deleted hand code to machine precision; RaviartThomas-type elements come out as pure +identity (Piola-equivalent) automatically. + +GuzmanNeilanFirstKindH1 (orders 0-2, 2D/3D) is also automatic: vertex/edge point +values are value point-groups mapping by $K$ (`_piola_point_rows`, the mirror of +`_scalar_point_rows`), and the trailing tangential facet constraints are dropped since +`PiolaBubbleElement.space_dimension()` already returns the reduced count (the Bell +pattern, inherited rather than passed as a parameter); `PiolaBubbleElement.__init__` +still provides the reduced `entity_dofs` bookkeeping. `GuzmanNeilanFirstKindH1` is +`class GuzmanNeilanFirstKindH1(PiolaPhysicallyMappedElement, PiolaBubbleElement)`: MRO +puts the automatic `basis_transformation` first, while `PiolaBubbleElement.__init__` +(reached via `super()`) still sets up `self._element` and the reduced dof bookkeeping. Its hand-derived vertex-facet coupling correction ("fix +discrepancy" in `finat/piola_mapped.py`) emerges automatically from the Vandermonde +residual elimination, since the per-point normal moments evaluate against vertex basis +functions of the extended element. + +**AlfeldSorokina and GuzmanNeilanH1div (divergence at vertices).** AlfeldSorokina's dual +basis interleaves, at every vertex, one `PointDivergence` node with `sd` `ComponentPointEvaluation` +nodes at the same point; GuzmanNeilanH1div wraps the same AlfeldSorokina dofs (plus GN +facet bubbles), so it needs the identical fix. `PointDivergence` does not fit the +existing `PhysicallyMappedFunctional` shapes: it has both `order = 1` (a derivative) and +`comp != ()` (a value component) on the *same* weight entries, in the pattern $\ell(v) = +\sum_i \partial_i v_i$ — a trace pairing an order-1 derivative multi-index with a rank-1 +value component, not a rank-one `direction` the SVD factorization can recover. Recovered +in `from_fiat` by testing, per point, that the $(alpha, comp)$ weight matrix is a scalar +multiple of the identity; represented with a new `divergence: bool` flag (order and rank +are otherwise meaningless for this node, so reusing them would have made `evaluate`/ +`pullback` silently wrong for it — a new field beats overloading old ones here). Handling +is a one-line closed form, not a generalization of the existing machinery: the +(contravariant) Piola pullback commutes with the divergence up to $\det J$, regardless of +which entity the node sits on, so `PiolaPhysicallyMappedElement._divergence_rows` strips +divergence nodes out of the group first and sets `V[i, i] = detJ * ell.weights[0]` +directly, before either `facet_dof_rows` or `point_dof_rows` runs. + +Stripping the divergence node still left a `facet_dof_rows` bug: AlfeldSorokina's edge +dofs (`dim == sd - 1` in 2D) are `ComponentPointEvaluation`s too — plain Cartesian point +values that happen to sit on a codimension-1 entity, not genuine facet moments. Routing +them through the normal/tangential frame logic built for MTW/JM/GN produced a singular +Gram matrix (`B @ B.T`), because that logic assumes the tangential-only residual group +has exactly `sd - 1` members (the facet completion count), but `sd` unrelated Cartesian +components don't reduce to that shape. Fix: `_is_cartesian_point_group` recognizes a +group of rank-1, order-0 nodes sharing one single point (real facet moments are built +from multi-point quadrature, even at low order — verified by checking `len(ell.points)` +for BernardiRaugel/MTW/JM/GN facet dofs, all $>1$) and dispatches it to the same +`_piola_point_rows` Cartesian-component transform `point_dof_rows` uses, regardless of +the entity's topological dimension. This is the same lesson as the FacetFrame work: +dispatch on what the node's *data* looks like, never on which entity it happens to sit on. + +Next steps: GN second kind (interior derivative moments need the divergence detJ rule, +same as above), ArnoldWinther (vertex tensor values + higher facet moments), the +extended-element path for reduced HCT (macro polynomial spaces); covariant elements. + +## Paper process (Kirby, Marsden & Brubeck, "Automation of Finite Element +Transformations"): staged plan + +Process agreed 2026-07-16 (full plan: `~/.claude/plans/lazy-yawning-quasar.md`): alternate +a *math stage* (derive/restate theory, hand-work an example, extract a glossary) with a +*code stage* (rename/refactor to match, verify via `test/finat/test_zany_mapping.py` + +`flake8` + `pydocstyle`), starting with the scalar case (Stages 0-3) and then testing the +hypothesis that Piola-mapped elements unify further by composing the Piola transform with +the scalar machinery (Stages 4-6), before assembling `fiat_zany_auto/paper.tex` (Stage 7). +Math work accumulates here before being promoted into the paper. Stage 0 (restate +pullback/push-forward duality, $M = V^T$, and $V = EV^cD$ precisely) is already satisfied +by the "Transformation Theory" section opening this file. + +### Stage 1 (2026-07-16): worked example — Hermite vertex jets + +The affine-interpolation-equivalent case (`ScalarPhysicallyMappedElement.point_dof_rows`, +`finat/zany.py:420-458`), worked by hand and checked to machine precision against +`finat.Hermite.basis_transformation` on a concrete triangle. + +**Setup.** Reference cell $\hat K$ = UFC triangle, vertex $\hat v_0 = (0,0)$. FIAT's cubic +Hermite dual basis carries, at $\hat v_0$, one value node and two derivative nodes +$\hat\ell_1 = \partial/\partial\hat x_1$, $\hat\ell_2 = \partial/\partial\hat x_2$ +(`FIAT.hermite.CubicHermite(ref).dual_basis()`, confirmed via `deriv_dict`). Physical +triangle vertices `((0.0, 0.1), (1.17, -0.09), (0.15, 1.84))` (the `phys_el` fixture of +`test/finat/conftest.py`), giving the constant Jacobian +$$J = \begin{pmatrix} 1.17 & 0.15 \\ -0.19 & 1.74\end{pmatrix}$$ +(`coordinate_mapping.jacobian_at`, i.e. $x = J\hat x + b$ — this is the *code's* $J$, +mapping reference to physical; see the reconciliation note below for how it relates to +the "paper's $J$" of the opening Transformation Theory section). + +**Mechanism.** Away from a facet there is no geometric frame (`FacetFrame` does not +apply), so `point_dof_rows` treats the whole order-1 group at $\hat v_0$ as its own +completion, exactly the affine-interpolation-equivalent case of Kirby (2017): the group +must span the derivative jet, and $V$'s $2\times2$ block on $\{1,2\}$ is obtained by +expanding each node's *pulled-back* direction in the group's own direction basis. +Precisely (`finat/zany.py:437-458`): + +1. `PhysicallyMappedFunctional.from_fiat` recovers each $\hat\ell_i$'s direction/weight + pair $(d_i, w_i)$ from the reference dual basis by a rank-1 SVD of the derivative + weight matrix — only determined up to a common sign/scale, so $w_i d_i$ (not $d_i$ + alone) is the invariant quantity, the *true* reference direction $e_i$. Here FIAT + stores $\hat\ell_1$ directly as $(d_1,w_1)=((1,0),\,1)$, but recovers $\hat\ell_2$ as + $(d_2,w_2) = ((0,-1),\,-1)$ (sign flipped by the SVD) — with $e_2 = w_2 d_2 = (0,1)$ + the correct Cartesian direction is recovered regardless. +2. `directions = [d_1; d_2]`, `Dinv = inv(directions.T)`: because $d_i$ (not $e_i$) enters + here, $Dinv$ absorbs the same sign ambiguity that $w_i$ will later correct. +3. For each pair $(i,j)$: `s = _weight_ratio(w_i, w_j)` $= w_i/w_j$ recovers the *relative* + sign/scale between the two recovered factorizations (`_weight_ratio` raises if the + weights are not simply proportional — here trivially $\pm1/\pm1$), and + `x = s * Dinv[col_j]`; then `V[i,j] = pullback(J, d_i) @ x`, where + `pullback(J, d_i) = J @ d_i` for an order-1 direction (verified numerically: the chain + rule for a first derivative is a plain linear map, no transpose, since `direction` is + contravariant in the *entity* Jacobian at this stage — the transpose appears only once + $d_i$ is re-expressed in the *dual* frame `Dinv` provides). + +Composing steps 2-3 for the true (sign-corrected) directions $e_i$ collapses to the clean +statement $V[i,j] = e_j^T J\, e_i$ — i.e. on this basis $V\big|_{\{1,2\}} = E^TJE$ with +$E = [e_1\,|\,e_2]$ (here $E=I$, the Cartesian basis, so $V\big|_{\{1,2\}} = J^T$ exactly); +the $w_i/d_i$ split and `_weight_ratio` correction exist only so this holds *regardless* +of which sign/scale `from_fiat`'s SVD happens to recover, not because the mathematics +needs a sign at all. + +**Verification.** Computed by hand from the formula above and checked bit-for-bit against +`finat.Hermite(ref).basis_transformation(mapping)` (evaluated via `gem.interpreter`): +$$ +V\big|_{\{1,2\}} = J^T = \begin{pmatrix} 1.17 & -0.19 \\ 0.15 & 1.74 \end{pmatrix} +\quad\Longleftrightarrow\quad +M\big|_{\{1,2\}} = V^T\big|_{\{1,2\}} = J = \begin{pmatrix} 1.17 & 0.15 \\ -0.19 & 1.74 \end{pmatrix}, +$$ +matching `M[1:3,1:3]` from the running code exactly (both diagonal *and* off-diagonal +entries, so the sign/scale bookkeeping above is exercised nontrivially, not just checked +on a symmetric or diagonal case). + +**Reconciliation with the opening Transformation Theory section.** Line 15 there states +the vertex-gradient block of $M$ as "$J^{-T}$ ... in the *paper's* $J$." The code's +`coordinate_mapping.jacobian_at` returns $J_{code}$ with $x = J_{code}\hat x + b$ +(reference $\to$ physical), whereas Kirby (2017)'s $F$ maps physical $\to$ reference, so +its Jacobian is $J_{paper} = J_{code}^{-1}$; hence $J_{paper}^{-T} = J_{code}^T$, exactly +the $M$ block computed above. **Glossary note for Stage 3**: every code docstring should +state explicitly which of these two (inverse, transpose-of-each-other) conventions a +given $J$ is, since the two papers' $F$ and the code's `coordinate_mapping` point in +opposite directions. + +### Stage 2 (2026-07-16): worked example — Morley edge completion + +The general (non-equivalent) completion case (`ScalarPhysicallyMappedElement.facet_dof_rows`, +`finat/zany.py:368-418`, using `FacetFrame`), worked by hand for edge 0 (the hypotenuse, +joining vertices 1 and 2) on the same physical triangle as Stage 1, and checked bit-for-bit +against `finat.Morley`. + +**Setup.** FIAT's Morley dual basis has, on edge 0, a single node $\hat\ell_3$: the +normal derivative at the midpoint $(0.5,0.5)$, stored as an *average* +(`deriv_dict` weights $(1,1)/\sqrt2$ on $(\partial_{\hat x_1}, \partial_{\hat x_2})$, +i.e. direction $\hat n = (1,1)/\sqrt2$, the reference normal itself — Morley has no +separate tangential-derivative node, so $\hat\ell_3$ must be *completed* using the +vertex value nodes ($\hat\ell_1,\hat\ell_2$ at $v_1=(1,0)$, $v_2=(0,1)$, already +processed since $\mathrm{order}=0$) before it can be expanded on the physical cell. +Reference edge-0 tangent (FIAT) $\hat t = (-1,1)$ (the chord $v_2-v_1$, unnormalized). + +**Mechanism.** `FacetFrame(Mo, 0, J)` builds the reference frame +$[\hat n\,|\,\hat t\,]$ and its physical image $[C\,|\,J\hat t\,]$ with +$C = $ `generalized_cross`$(J\hat t)$ (here just a $90°$ rotation of $J\hat t$ in 2D). +`facet_dof_rows` then, for $\hat\ell_3$: + +1. `reference_coefficients(\hat\ell_3.direction)` solves $\hat d = a\hat n + \beta\hat t$ + for $(a,\beta)$ in the *reference* frame — here $a=-1$ (sign from the SVD recovery + of `from_fiat`, harmless: every downstream quantity built from $a$ is consistently + rescaled), $\beta=0$ exactly, since $\hat\ell_3$'s direction *is* the reference normal + with no tangential component by construction. +2. `decompose(\hat\ell_3.pullback(J).direction)` solves the *physical*, symbolic system + $J\hat d = x_0 C + x_1 J\hat t$ via `adjugate`/`determinant` (Cramer's rule, since $J$ + is symbolic in general; numeric here because $J$ is the constant matrix of Stage 1) + giving $(x_0,x_1) = (0.612629,\,-0.244111)$. +3. $c = x_0\cdot(\|C\|/\kappa)/a = 1.337343$ becomes $V[3,3]$ — the diagonal + "own-node" coefficient, i.e. the $B_{nn}$ entry of the notebook's $2\times2$ block + $B = \hat G J^{-T} G^T$ (line 34 of the Transformation Theory section above; that + block's off-diagonal $B_{nt}$ entry is exactly $-c\beta/a$-type bookkeeping, here + $0$ since $\beta=0$). +4. The tangential residual $r = x_1 - c\beta = -0.244111$ is *not itself* a new + unknown: because $\hat t$ is a **mapped reference tangent**, the reference + functional "derivative along $\hat t$ at the midpoint" coincides with a functional + already expressible in the element's own reference basis — `ell.with_direction(\hat t) + .evaluate(Mo)` (a numeric generalized-Vandermonde row) gives coefficients + $(0,1,-1,0,0,0)$: the classical fact that for a quadratic, the derivative along the + full chord at the midpoint equals $f(v_2)-f(v_1)$ exactly (Brubeck & Kirby 2025's + univariate-exactness/FTC argument realizing $D$, here in its simplest 1-node form). +5. Elimination: $V[3,1] \mathrel{+}= r\cdot1 = -0.244111$, $V[3,2] \mathrel{+}= r\cdot(-1) + = +0.244111$ — both already-`processed` vertex rows, so no `NotImplementedError`. + +**Verification.** `finat.Morley(ref).basis_transformation(mapping)` (same physical +triangle as Stage 1) gives column 3 of $M=V^T$ as +$M[:,3] = (0,\,-0.244111,\,0.244111,\,1.337343,\,0,\,0)$, matching $V[3,\cdot]$ from +steps 3-5 to 6 decimal places. + +**What this example newly exercises, vs. Stage 1.** Stage 1's Cartesian point-jet +group needed only $E=I$ (no frame) and a diagonal-plus-off-diagonal $V^c$ block coming +purely from the chain rule. Here, for the first time, $E$ (the frame's normal/tangential +extraction, keeping only the normal row) and $D$ (the numeric elimination of the +tangential residual through *already-assembled, lower-dimension* rows — vertex nodes, +processed before edges by the increasing-dimension loop of `basis_transformation`) are +both genuinely nontrivial and distinct, which is exactly the completion mechanism the +notebook's $V=EV^cD$ factorization (lines 21-36) describes in the abstract. + +The implementation mirrors the theory factor by factor: + +1. **Sparsity from topology**: each row of $V$ is a reference node; its nonzero columns + are the physical nodes on the same entity (via $V^c$) plus the nodes on adjacent + entities pulled in by $D$ (e.g. edge rows couple to the edge's vertex dofs, and for + hierarchical elements to the same edge's trace moments). Derive this from + `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. +3. **$D$ from unisolvence, not from closed-form rules**: the univariate exactness rules + (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. +4. Constrained spaces (reduced HCT, Bell) reuse the same machinery on the extended + element; the FIAT element must expose the constraint functionals as extra dofs. + +### Stage 3 (2026-07-16): scalar glossary + code refactor — done + +Deliverables, applied to `finat/zany.py` only (zero semantic diff, no behavior change): + +* A module-level glossary docstring mapping every recurring symbol (`V`, `E`, `V^c`, `D`, + `J`, `K`, `normal`, `direction`/`weights`/`order`/`rank`/`divergence`, `a`, `beta`, `x`, + `c`, `r`, `Dinv`, `B`/`Binv`/`L`/`Lmap`) to its code location and mathematical role, and + stating explicitly which Jacobian-inverse convention `J` uses relative to Kirby (2017)'s + $F$ (the Stage 1 finding). +* `_materialize_jacobian(J)`: a new free function replacing three copies of the same + `numpy.array([[J[i, k] ...` snippet (`PiolaFacetFrame.__init__`, `_divergence_rows`, + `_piola_point_rows`), so the numpy-materialized Jacobian has one name and one + construction site instead of three silently-identical ones. +* `FacetFrame.normal`/`PiolaFacetFrame.normal`'s differing scaling convention + (`compute_normal` vs `compute_scaled_normal`) is now stated in the module glossary, + not just each class's own docstring, so it's visible without reading both classes. + +Verified: `pytest test/finat/test_zany_mapping.py test/finat/test_zany_automation.py` +(124 passed), the full `test/finat/` suite (338 passed, 8 skipped), `flake8 finat/zany.py`, +`pydocstyle finat/zany.py` — all clean, no matrix entries changed. + +This closes the first math-code loop (Stages 0-3, scalar case). Stages 4+ (Piola +composition hypothesis and beyond) remain only roughly sketched in the plan file and will +be planned in detail next. + +### Paper migration (2026-07-16) + +The plan now lives at `~/git/fiat_zany_auto/PLAN.md`, and the scalar theory +(Stages 0-3 above) has been drafted as Section 2 of `~/git/fiat_zany_auto/paper.tex` +("Automation of finite element transformations", Kirby, Marsden & Brubeck; acmart, +compiles cleanly with `make`). Key notation choices made in the draft, to keep +consistent going forward: $F: K \to \hat K$ physical-to-reference (Kirby 2017), but $J$ +defined by $F^{-1}(\hat x) = J\hat x + b$ (the *code's* Jacobian, stated explicitly in +eq. 2.3); the generic node $\ell_{X,w,d}$; the push-forward closed form +$F_*(\ell_{X,w,d}) = \hat\ell_{\hat X, w, (J^{-1})^{\otimes m}d}$; the "design equation" +$J^{\otimes m}\hat d = \sum_j c_j d_j$ as the unifying statement of every row solve; and +a theory-to-code correspondence table (Table 1). The Piola/implementation/examples +sections are `\todo` stubs pending Stages 4+. + +### Review response, RCK round 1 (2026-07-17) + +Robert's remarks (`rck_remarks_for_claude.md`) addressed in `~/git/fiat_zany_auto/paper.tex` +(now 12 pp., compiles clean). Disposition: + +- **det J > 0 dropped everywhere.** §2.1 now states both orientations occur under UFC; + the normal formula (2.14) is presented as orientation-covariant (κ = ±1, ν flips with + C exactly as the reference formula would on a reflected cell); (2.23) notes Gram + determinants are positive and det J enters linearly. *Verified in code*: ran + `check_zany_mapping` on a clockwise triangle (det J = −2.06) for Morley, Hermite, + Bell, Argyris(5,6,7, avg=True), WuXuH3NC, WuXuRobustH3NC, and on a reflected tet + (det J = −1.28) for Morley/Hermite — all pass. **Note: the test-suite fixtures only + use positively oriented cells; suggest adding a flipped-cell parametrization.** +- **Notation**: physical/reference normals renamed N, n̂ → ν, ν̂ (avoids clash with + nodes n̂_j and node-set 𝒩). Facet frame and scaled tangents (t̂_k = v̂_k − v̂_0) + now defined explicitly (new eqs. framedef/scaledtangents). +- **(2.8) scope**: one direction + one order per functional stated as deliberate, + no premature generalization; symmetry-WLOG + rank-deficiency remark (d = ν⊗ν rank 1; + construction never uses rank; covers Wu–Xu-type nodes). +- **"dividing by J"** → multiplying each slot by J^{-1} (gradients by J^{-T}); + invariance/design-equation paragraphs rewritten ("conversely" contrast made explicit; + numeric-vs-symbolic split cued early); no-type-dispatch point added ("a normal + component is detected through a ≠ 0, not declared through a class hierarchy"). +- **Binning**: §2.4 now states grouping = entity (element data) + order (numeric); + shared points/weight-ratios/direction-basis are checked numerically, nothing else needed. +- **FTOC proof**: new Prop. 2.1 (exactness identities are recovered — unisolvence + uniqueness, 1-line proof), with FTOC instance proved (average → all C¹; Morley + midpoint → on P̂ via linearity of ∂_t f along the edge) and IBP/Stokes-type (Wu–Xu + face moments) as further instances. +- **E, V^c, D separated then recursion (his 2-part request)**: §2.6 now defines the + three factors exactly (V^c block per eq:Vcblock, E Boolean, D rows = w·V — self- + referential but triangular by entity structure), §2.7 "Row recursion" evaluates the + product without materializing (graph/topological-order remark included), with + **Algorithm 1** (algpseudocode) mirroring `basis_transformation`/`facet_dof_rows`. +- **Worked examples**: explicit 6×6 Morley V with closed forms c_e = det J |t̂|/|Jt̂|, + r_e = (Jν̂·Jt̂)/|Jt̂|² — *verified to 6.7e-16 against the code on both orientations*; + hexic Argyris completion rows computed: ℓ̂⁽⁰⁾ = δ_{v+} − δ_{v−} (FTOC) and + ℓ̂⁽¹⁾ = 3δ_{v+} + 3δ_{v−} − μ̂_ê (IBP, couples to same-edge trace moment; exact + integers to machine precision). +- **Novel-element example**: sec:examples todo now points at WuXuH3NC/WuXuRobustH3NC, + which the automation already handles (tested above). +- **Honesty fix**: scalar completions only ever couple to *identity* rows; the fully + recursive regime (coupling into non-identity rows) is now correctly attributed to the + Piola macroelements (§3, deferred). +- New bib entries: `ufc` (Alnæs–Logg–Mardal FEniCS-book chapter), `wu-xu` (Math. Comp. 2019). + +Not done (open for co-authors): whether more calculations deserve proposition status +(only Prop 2.1 added so far); explicit Stokes-identity theorem for tet face moments +(covered as an instance of Prop 2.1, no standalone theorem); India's email still TODO. + +### Stage 4 (2026-07-17): Piola theory derived as a composition with the scalar machinery + +Goal (user hypothesis from PLAN.md): derive the Piola-mapped case *from* the scalar +case by composing basis transformations, as the route to unifying +`ScalarPhysicallyMappedElement` and `PiolaPhysicallyMappedElement`. Result: the +hypothesis is correct, the composition is exact, and it was verified to machine +precision against every Piola element in the FInAT test zoo (automated *and* +hand-coded), in 2D and 3D, on both orientations. Details below; conventions as in +the paper (`~/git/fiat_zany_auto/paper.tex` §2): $F : K \to \hat K$, +$F^{-1}(\hat x) = J\hat x + b$. + +**1. Composition theorem.** On an affine cell every Piola pullback factors through +the componentwise scalar pullback $F^*_0(\hat f) = \hat f \circ F$: +$$F^*_{\mathrm{piola}} = \theta_A \circ F^*_0, \qquad (\theta_A f)(x) = A\,f(x),$$ +with $A$ a *constant* matrix acting on the value components: $A = J/\det J$ +(contravariant), $A = J^{-T}$ (covariant), slot-wise for double variants. +Dually, on functionals, +$$F_*^{\mathrm{piola}} = F_*^0 \circ \theta_A^*, \qquad \theta_A^*(n) = n \circ \theta_A,$$ +and both factors act *slot-locally* on the generic node: extend the scalar node to +$\ell_{X,W}(f) = \sum_q \langle W_q, \nabla^m f(x_q)\rangle$ where $W_q$ now carries +$r$ value slots and $m$ derivative slots. Then the unified push-forward is +$$F_*(\ell_{X,W}) = \hat\ell_{\hat X, \hat W}, \qquad +\hat W = (\text{each contravariant value slot}) \cdot \Theta^T + \;\otimes\; (\text{each covariant value slot}) \cdot J^{-1} + \;\otimes\; (\text{each derivative slot}) \cdot J^{-1},$$ +with $\Theta^T = J^T/\det J = K^{-1}$, $K = \det J \, J^{-T}$ the cofactor matrix. +Points and scalar weights are preserved (paper eq. 2.9 is the $r=0$ case). Note the +collapse: **covariant value slots transform exactly like derivative slots** — a +covariant element is "scalar theory with one more $J^{-1}$ slot". + +**2. Mirror duality of invariance.** A physical slot direction $\delta$ is invariant +iff $\delta = \rho^{-1}\hat\delta$ for the slot map $\rho$: +- derivative / covariant slots: $\delta = J\hat\delta$ — mapped tangents (scalar case, + paper §2.2); +- contravariant slots: $\delta = K\hat\delta$ — cofactor images. + +The geometric content is one lemma used by both cases: the generalized cross +product intertwines $J$ and $K$, $\bigwedge_k J\hat t_k = K \bigwedge_k \hat t_k$ +(i.e. $C = K\hat C$). In the scalar case this gave the physical unit normal +$\nu = \kappa C/|C|$; in the Piola case it gives the **cofactor lemma**: FIAT's +`compute_scaled_normal` is exactly $K$-covariant, +$$K \hat\nu^s = \nu^s_{\mathrm{phys}} \quad \text{(verified: error} \le 3\cdot10^{-16}, +\text{all facets, 2D+3D, both orientations)},$$ +so physical normal-flux moments are invariant *with no symbolic factor at all* +(the mirror of the scalar case's invariant mapped-tangential derivatives — but +exact, no $\kappa$ needed). Consequence: FIAT's physical scaled normal is NOT +always outward on negatively oriented cells; its sign rides the UFC vertex +ordering, which is what makes the invariance orientation-robust. + +**3. Facet frame, mirrored, with closed-form net.** Reference frame +$[\hat\nu^s \,|\, \hat t_k]$ (scaled normal + scaled tangents). FIAT's physical +facet-dof conventions (read off `FIAT/mardal_tai_winther.py` and verified across +the zoo): +- normal profiles: against $\nu^s = K\hat\nu^s$ → invariant; +- tangential profiles, 2D: against plain mapped tangents $J\hat t$; +- tangential profiles, 3D: against $\nu^s \times (J\hat b)$, cross products of the + scaled normal with mapped in-plane (RT-profile) vectors — *same formula on + reference and physical cells*, the same principle as the scalar normal. + +Pushing the 3D convention forward with the slot calculus, using the identity +$\Theta^T(Ka \times Jb) = (J^{-1}Ka)\times b$ and BAC-CAB, the net map on frame +coordinates is +$$N = \begin{pmatrix} 1 & -\det J\, \frac{\hat t_l \cdot M^{-1}\hat\nu^s}{|\hat\nu^s|^2} \\ 0 & s\,I_{d-1} \end{pmatrix}, +\qquad s = \det J\, \frac{\hat\nu^s \cdot M^{-1} \hat\nu^s}{|\hat\nu^s|^2}, +\qquad M = J^T J .$$ +The tangential block is a **scalar multiple of the identity** — the +reciprocal-basis correction $S^{-1}$ and the mixing matrix $Y$ of +`PiolaFacetFrame` are derived consequences of this, not independent structure; +no in-plane matrix correction is needed. The top-right entries are the +completion residuals (the mirror of the scalar $r_k$): a pulled-back tangential +profile leaves behind a *normal*-direction functional, eliminated numerically +through the generalized Vandermonde row (`evaluate` on the nodal basis), exactly +as the scalar case eliminates *tangential* residuals. Spot-verified numbers +(fixture cells): MTW 3D face 0: $s = 1.730944$, completion row +$(0.006240, -0.020316)$, matching an unconstrained least-squares fit of FIAT's +actual physical duals to 9e-16; GN2 2D edge 0 tangential dof: pushed frame +coords $(0.398465, 1.154215)$ = closed form. + +**4. Divergence nodes.** $\mathrm{div}$ = contraction of one contravariant value +slot with one derivative slot; the slot maps contract to +$(K^{-1})(J^{-1})^T \delta = \delta/\det J$, so +$F_*(\ell_{\mathrm{div}}) = (\det J)^{-1} \hat\ell_{\mathrm{div}}$ — invariant up +to the scalar $\det J$, independent of entity ("the divergence miracle", = +`_divergence_rows`). Tensor case (moments of $\mathrm{div}\,\sigma$ against +vector weights, the Arnold-Winther constraint nodes): one value slot survives +and maps by $\Theta^T$, same $1/\det J$; covered by the calculus, currently +blocked only by `from_fiat` not parsing `IntegralMomentOfTensorDivergence`. + +**5. Point data vs moments.** Single-point order-0 nodes (vertex values, edge +midpoint values, interior point values — C0-type data) keep Cartesian components; +their push-forward is the plain $\Theta^T$ slot map and the group block-solve +gives $K$-blocks per slot (the exact mirror of the scalar vertex jets, whose +blocks are $J^T$ per slot). Interior *moments* are invariant by convention +(physical tests are Piola-mapped). **Bug-relevant discovery**: the current +`PiolaPhysicallyMappedElement.invariant_dofs` rule (`order == 0 and dim == sd`) +misclassifies interior Cartesian *point* data as invariant — it happens to be +correct for every currently automated element, but is wrong for Guzman-Neilan +second kind (interior barycenter values → need $K$ blocks, as its hand-coded +transformation does) and Hu-Zhang point variant (interior point values of +$\sigma$). Discriminator that works for the whole zoo: single-point order-0 +nodes are point data; multi-point are moments. (HZ-point facet dofs additionally +show rank-2 single-point *frame* dofs, so the facet-level Cartesian test +`len(points)==1 and rank==1` is still needed there.) + +**6. Unified formulation = duality.** Everything above collapses into one +statement: with slot-mapped weights, the generalized Vandermonde matrix +$$B_{ij} = F_*(n_i)(\hat\psi_j)$$ +is computable for any pullback that factors through $\theta_A$, and $V = B^{-1}$, +inverted blockwise along the entity hierarchy — the row recursion of paper +§2.7 is precisely this block-triangular inversion. A ~60-line prototype +(numeric $J$; slot maps + frame net + divergence + point-data rules, nothing +else) reproduces `basis_transformation` to machine precision for: +MTW(1,2), JohnsonMercier, ArnoldWintherNC, AlfeldSorokina, BernardiRaugel(+Bubble), +GuzmanNeilan first kind (1,2), **GuzmanNeilan second kind**, GN Bubble, GN H1div, +ReducedArnoldQin, ChristiansenHu, **HuZhang 3/4 in both point and integral +variants** — 2D and 3D where defined, positive and negative orientation +(bold = elements the current automated class cannot handle; their hand-coded +transformations agree with the composition theory). All `check_zany_mapping` +ground-truth tests also pass on reflected cells (fixtures only test positive +orientation; worth adding). + +**7. What this means for Stage 5 (code unification).** The +`{Scalar|Piola}PhysicallyMappedElement` split reduces to a small table keyed by +the FIAT mapping string: + +| ingredient | affine (scalar) | contravariant piola | covariant piola | +|-----------------------------|-----------------------------|------------------------------|-----------------| +| value slot map | $I$ | $\Theta^T = K^{-1}$ | $J^{-1}$ | +| derivative slot map | $J^{-1}$ | $J^{-1}$ | $J^{-1}$ | +| frame-stable directions | tangents ($J\hat t$) | scaled normal ($K\hat\nu^s$) | tangents | +| completed directions | unit normal $\nu$ | in-plane (2D: $J\hat t$; 3D: $\nu^s{\times}J\hat b$) | normal | +| point-data block | $J^T$ per slot | $K$ per slot | $J^T$ per slot | +| special contractions | — | div → $\det J$ | curl analog | + +plus one shared engine: recover node data numerically (`from_fiat` extended with +value slots — already done), decide invariance from frame profiles, expand the +non-stable part symbolically, eliminate residuals through `evaluate`, recurse +over entities. `FacetFrame`/`PiolaFacetFrame` merge into one class parametrized +by the stable subspace; the $S^{-1}$/`Y` machinery is replaced by the closed-form +$N$ above. Newly automatable for free: GN2, HuZhang (both variants), BR, CH, RAQ, +and AW/AWnc once `from_fiat` learns tensor-divergence moments. Covariant (Nedelec +-type zany) elements: theory ready, value slot $\equiv$ derivative slot. + +Verification artifacts: the prototype now lives in the repo as +`test/finat/test_claude_zany_piola.py`, upgraded from the in-session scripts: +it computes $V = B^{-1}$ by *block back-substitution* (no dense inversion +anywhere) and asserts both the support law of $B$ (see the Stage 5 design +below) and agreement with `basis_transformation`, for the whole zoo, 2D and +3D, both orientations (56 parametrized cases). + +### Stage 5 design (2026-07-17): sparsity, cost, and the unified elimination + +Robert's (correct) objection to the first prototype: it formed $B$ densely and +called `inv`. That was a verification shortcut, not the algorithm. This section +records the sparsity theory that makes the duality formulation *no less sparse +and no more expensive* than the current hand-structured code, and the +implementation design in enough detail to be reviewed before touching +`finat/zany.py`. + +**1. Support law (sparsity theorem for $B$).** Order the dofs by the entity +partial order: $(d, e) \preceq (d', e')$ iff $e \subseteq +\overline{e'}$ (closure). Claim: $B_{ij} = F_*(n_i)(\hat\psi_j) \ne 0$ only if +$\mathrm{entity}(j) \preceq \mathrm{entity}(i)$, so $B$ is block lower +triangular with per-entity diagonal blocks. Sharper, by dof kind: + +* *Single-point rows* (vertex/edge/interior Cartesian point data): exactly + block-diagonal within the point group. Proof: $F_*(n)$ at a point $\hat x$ is + the $\Theta^T$-combination of the component-evaluation nodes at $\hat x$ — a + pointwise identity of functionals on $C^0$, needing no polynomial exactness. +* *Divergence rows*: $B_i = e_i/\det J$ exactly. Proof: the divergence miracle + $(K^{-1}\!\otimes\!J^{-1})\delta = \delta/\det J$ turns $F_*(\ell_{div})$ + into $\hat\ell_{div}/\det J$, and $\hat\ell_{div}$ *is* a reference node, so + its generalized Vandermonde row is $e_i$. +* *Interior invariant moments*: $B_i = e_i$ by definition of invariance. +* *Facet rows*: support $\subseteq$ facet block $\cup$ dofs on the strict + closure $\partial f$. Proof sketch: decompose the slot-mapped weight in the + reference frame; the frame components give the facet-block entries (the net + $N$); the residual $R$ is a moment supported on $f$, hence a functional of + the trace of its argument on $f$. For $j \notin \overline{f}$, $\hat\psi_j$ + has vanishing dofs on $\overline{f}$, and *trace unisolvence* — the dofs on + $\overline{f}$ determine the traces on $f$ that facet functionals sense, + i.e. exactly the property that makes the element (H1- or H(div)-) + assemblable across facets — forces $R(\hat\psi_j) = 0$. This is the Piola + analogue of the scalar exactness identities (paper Prop 2.1): in Morley, + "the trace on $e$ is determined by closure dofs" is implemented by FTOC + along the edge. **Sparsity of $B$ = conformity of the dof layout.** + +Verified numerically (scratch script, 2026-07-17): max $|B_{ij}|$ outside the +predicted support is $0$ to machine precision for all 28 element instances of +the zoo, on a negatively oriented cell; the test file asserts the same law on +both orientations on every run. + +**2. Elimination = block back-substitution = the paper's row recursion.** +$BV = I$ with $B$ block triangular gives, per entity $e$ with closure dof set +$c$ (already processed, since dims are visited in increasing order): +$$V_e = B_{ee}^{-1}\,(I_e - B_{ec} V_c).$$ +This is the scalar row recursion $V_i = c\,e_i + \sum_{j \in P} \gamma_{ij} +V_j$ in matrix clothing. Fill-in: $\mathrm{supp}(V_e) \subseteq e \cup c$, +because the closure relation is *transitive* — the predicted fill equals the +final fill, nothing grows during the elimination. Depth of the recursion +$\le sd$ (vertex → edge → face → interior). At no point is a dense matrix +formed or inverted; the prototype's `np.linalg.inv` is gone +(`composition_transformation` in `test/finat/test_claude_zany_piola.py` now +implements exactly this loop and asserts the support law en route). + +**3. Diagonal blocks invert in closed form (the symbolic cost).** The only +$J$-dependent inverses ever needed are: + +* point group: inverse slot map $= K$ per slot (already what + `_piola_point_rows` emits; $K = \mathrm{adj}(J)$ is *polynomial* in $J$); +* divergence: $\det J$ (polynomial); +* facet frame: $N$ is block triangular + $\begin{pmatrix}1 & r_l\\ 0 & sI\end{pmatrix}$ with **scalar** diagonal, so + $N^{-1} = \begin{pmatrix}1 & -r_l/s\\ 0 & s^{-1}I\end{pmatrix}$, and the + coefficients are Gram contractions of the $K$-mapped frame — no symbolic + $M^{-1} = (J^TJ)^{-1}$ solve, since $J^{-T}\hat\nu = K\hat\nu/\det J$: + $$s = \frac{|K\hat\nu^s|^2}{\det J\,|\hat\nu^s|^2} + = \frac{|\nu^s|^2}{\det J\,|\hat\nu^s|^2}, \qquad + r_l = -\frac{(K\hat t_l)\cdot(K\hat\nu^s)}{\det J\,|\hat\nu^s|^2}$$ + (verified to 1e-15, both orientations, 2D and 3D). Note the mirror of the + scalar edge coefficient $c_e = \det J\,|\hat t_e|/|J\hat t_e|$: numerator + and denominator swap roles under $J \leftrightarrow K$, tangent + $\leftrightarrow$ normal — the covariant/contravariant duality again. + +Everything else in $B_{ee}$ and $B_{ec}$ is *numeric and $J$-independent*: +restricted generalized Vandermonde data (reference tabulations contracted +with reference weights), computable and cacheable at element construction. +Each symbolic entry of $V_e$ is a sum of $\le$ (number of frame monomials) +terms (numeric constant) $\times$ (monomial in $\{1, s, r_l, \det J, K_{ab}\}$) +— for rank-$r$ elements the monomials come from $N^{\otimes r}$, still a fixed +finite set. This is the paper's $V = E\,V^c D$ factorization surfacing again: +numeric completion $E$, symbolic (block-)diagonal $V^c$, diagonal scaling $D$. +Per-dof symbolic cost: $O(sd^2)$ rational scalars; total GEM size $O(n\,sd^2)$ +— no symbolic LU, and *cheaper* than the current `PiolaFacetFrame` path, which +builds a symbolic physical Gram $G$ and its determinant for `Sinv`. + +**4. Implementation plan (what actually changes, reviewable pieces).** + +1. `finat/functional.py`: give `PhysicallyMappedFunctional` typed value slots + (mapping string per slot: affine / contravariant / covariant) and one + `pushforward(J, K)` that applies $J^{-1}$ per derivative slot (the existing + `pullback`) and $\Theta^T$ resp. $J^{-1}$ per value slot. Extend `from_fiat` + to tensor-divergence moments (trace over one derivative+value slot pair, + remaining value slot kept) — unlocks AW/AWnc. +2. `finat/zany.py`: one `basis_transformation` loop over entity blocks in + increasing dimension (the skeleton already exists); per-block behavior + chosen from *data* (divergence flag, single-point, entity dim, mapping + string), not from the Scalar/Piola subclass. `FacetFrame` and + `PiolaFacetFrame` merge into one frame class parametrized by the stable + subspace; it precomputes the numeric restricted Vandermonde/profile data at + construction and emits the $\{c\ \text{or}\ s, r_l, \det J\}$ closed forms + at transformation time. Delete `Y`/`Sinv`; `_piola_point_rows` and + `_divergence_rows` become instances of the generic slot-map inverse. +3. Fix the interior point-data rule (`invariant_dofs`): single-point order-0 + nodes are point data (K blocks), multi-point are invariant moments — + unlocks GN2 and HuZhang-point. +4. Runtime guard = the support law: assert the numeric residual expansion + vanishes outside the closure block (reference-time, once per element), so + a non-conforming dof layout fails loudly instead of producing a dense or + wrong $V$ (the generalization of Algorithm 1's assert line). +5. Tests: keep `test_claude_zany_piola.py` as the independent numeric ground + truth (it shares no code path with `basis_transformation`); add + reflected-cell fixtures to `test_zany_mapping.py`. + +Risk framing unchanged from PLAN.md: abandon any piece that adds more +complexity than it removes; the theory above is the review artifact to be +signed off *before* production code moves. + +--- + +## Stage 5 SHIPPED — duality engine as the default `basis_transformation` (2026-07-19) + +The mixin-unification plan above (Stage 5c/"Implementation plan" items 1–3) +was **superseded by a course change**: instead of merging the frame classes, +the numeric prototype of `test/finat/test_claude_zany_piola.py` was turned +directly into product code. `PhysicallyMappedElement.basis_transformation` +(`finat/physically_mapped.py`) now assembles $B_{ij} = n_i(\hat\psi_j \circ +F^{-1})$ symbolically in GEM and inverts it sparsely; `finat/zany.py` shrank +from ~900 lines (FacetFrame, PiolaFacetFrame, four hooks, two mixins) to two +thin configuration mixins (`_check_mapping` guards + the Piola `dof_scale` +override). Paper §4 ("Assembly by duality and sparsity-preserving +inversion") documents the theory. + +**Engine design (all in `finat/physically_mapped.py`):** + +1. *Per-slot maps, adjoint on the tabulation.* One effective matrix per + slot: $G = J^{-1}\Phi$ (derivative slots), $G = (J/\det J)^T\Phi$ (value + slots), applied to the reference direction/weights via the generalized + `PhysicallyMappedFunctional.pullback(A)` (now accepts any numeric or GEM + matrix); rows of $B$ are then numeric Vandermonde pairings. Numerator/ + denominator form keeps numerators polynomial: $s_i = (\det J\,|K\hat n|)^m$ + for scalar facet rows, $\det J^{2r}$ for Piola facet rows (using the exact + identity $J^T K\hat\nu^s = \det J\, \hat\nu^s$), $\det J$ for divergence. + +2. *Sparsity and constants inferred from numerics.* $B$ is assembled at 3 + deterministic sample Jacobians (`_sample_jacobians`, cosine-based, one + orientation-reversing). Entries below `pattern_tol = 1e-10` (row-relative) + at every sample are structural `gem.Zero`; entries equal at every sample + are constants (snapped to integers where they round), so invariant rows + never touch symbolic algebra. No closure bookkeeping, no invariance case + analysis, no support-law *assumption* — the pattern is measured. + +3. *SCC block triangularization + sparse back-substitution.* Tarjan on the + coupling digraph emits diagonal blocks in dependency order; + $V_S = (B'_{SS})^{-1}(\mathrm{diag}(s_i) I_S - B'_{Sc}V_c)$ with one + adjugate/determinant pair per block, cached across identical blocks + (vertex jets of different vertices hit the cache). GEM constant folding + through `Zero` keeps everything sparse; gotcha: numpy object arrays must + be the *left* operand of elementwise products, else GEM broadcasts a + vector node (`V[j] * coef`, not `coef * V[j]`). + +**Verification:** full finat suite 430 passed / 8 skipped (includes +`test_zany_mapping` direct physical fits, both orientations, and the +independent prototype cross-check in `test_claude_zany_piola.py`, which +still shares no code path with `basis_transformation`); flake8 + pydocstyle +clean; opaque `gem.Variable` Jacobian path verified bit-identical to the +literal path with max 3–9 nonzeros per row of $M$. + +**Numeric groundwork:** pattern-only elimination validated over the entire +scalar + Piola zoos before the symbolic build (scratch `validate_pattern.py`): +max diagonal block 5×5 (BrambleZlamalC2 order-4 jets), pattern leak ≤ 1e-14 +except BZC2 (~1e-10, its known FIAT-side conditioning). Hand-coded +transformations (HCT, PS, C2 elements, WuXu, Walkington) still override the +engine and are untouched.