Skip to content

Give each pull back an inverse - #71

Closed
pbrubeck wants to merge 3 commits into
pbrubeck/interpolate-holesfrom
pbrubeck/pullback-apply-inverse
Closed

Give each pull back an inverse#71
pbrubeck wants to merge 3 commits into
pbrubeck/interpolate-holesfrom
pbrubeck/pullback-apply-inverse

Conversation

@pbrubeck

@pbrubeck pbrubeck commented Sep 8, 2026

Copy link
Copy Markdown

Description

A pull back maps a function on the reference cell to the physical cell. UFL had no map in the other direction, so a form compiler that has to evaluate a physical expression on the reference cell — to dual-evaluate an Interpolate, for instance — carried its own copy of the Piola maps. TSFC has one.

This adds apply_inverse to AbstractPullback and implements it for every pull back that has an inverse:

pull back apply apply_inverse
IdentityPullback g g
ContravariantPiola (1/detJ) J g detJ K g
CovariantPiola Kᵀ g Jᵀ g
L2Piola g / detJ g detJ
DoubleContravariantPiola (1/detJ)² J g Jᵀ detJ² K g Kᵀ
DoubleCovariantPiola Kᵀ g K Jᵀ g J
CovariantContravariantPiola (1/detJ) Kᵀ g Jᵀ detJ Jᵀ g Kᵀ
MixedPullback per sub-element per sub-element
SymmetricPullback per unique block per unique block

PhysicalPullback and CustomPullback are their own inverse. UndefinedPullback inherits the NonStandardPullbackException from AbstractPullback, as it does for apply.

The test asserts the property that matters rather than the formulas: pushing a reference function forward and mapping it back returns the original, once remove_component_tensors and cancel_jacobian_products have cancelled the Jacobian products the round trip leaves behind. That closes exactly for eight of the nine, mixed and symmetric elements included. The L2Piola round trip leaves detJ / detJ standing, because cancelling a scalar factor is not something cancel_jacobian_products does, so that one is asserted directly instead.

Stacking

pbrubeck/pullback-apply-inverse sits on pbrubeck/interpolate-holes (FEniCS#497), so that it pairs with pbrubeck/form-interp-tsfc in Firedrake. That one branch is also the head of FEniCS#511, which proposes the same commit upstream; FEniCS#497 is in that diff for the same reason.

The pull back work itself does not depend on FEniCS#497 and can be read on its own.

Nothing calls apply_inverse yet. #69 uses it to lower interpolation in compute_form_data, which lets TSFC delete its copy of the maps.

🤖 Generated with Claude Code

https://claude.ai/code/session_0155wQMTMSGiF9a6TbyTVK5P

@pbrubeck
pbrubeck force-pushed the pbrubeck/interpolate-holes branch from 4f1550b to a75e36a Compare September 8, 2026 10:51
@pbrubeck
pbrubeck force-pushed the pbrubeck/pullback-apply-inverse branch 2 times, most recently from 6f6b47d to 3a43502 Compare September 10, 2026 23:09
A pull back maps a function on the reference cell to the physical cell.
Nothing in UFL went the other way, so a form compiler that has to evaluate a
physical expression on the reference cell had to carry its own copy of the
maps.

apply_inverse is the inverse of apply, for every pull back that has one.
PhysicalPullback and CustomPullback are their own inverse, and
UndefinedPullback inherits the refusal from AbstractPullback.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0155wQMTMSGiF9a6TbyTVK5P
@pbrubeck
pbrubeck force-pushed the pbrubeck/pullback-apply-inverse branch from 3a43502 to fbdb7d7 Compare September 10, 2026 23:40
pbrubeck and others added 2 commits September 11, 2026 08:13
* Fix Interpolate holes

* Interpolate._cache

* Expose Interpolate to the form compiler

* Fix shape and negation

* ruff

* Fix mypy

* Give a base form operator all of its domains

A base form operator's argument slots are not among its operands, so
traversing the operands misses the domains they are defined over. An
interpolation onto a point cloud is the case that notices: the target is
reachable only through the dual argument, and the form compiler numbers
its domains against this list.

extract_unique_domain still answers with the operand's domain, which is
where the operator takes its value.

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

* ruff

* Give every base form operator its domain analysis

BaseForm requires _analyze_domains and ufl_domains of its subclasses, and
BaseFormOperator never supplied them: Interpolate carried the pair alone, so
every other operator fell through to the deprecated Expr.ufl_domains, which
answers by calling extract_domains and recurses forever now that
extract_domains routes a base form operator back through it. Differentiating
an interpolation is the case that notices, as BaseFormOperatorDerivative has
no such override.

The analysis reads ufl_operands and argument_slots, both of which every base
form operator has, so it belongs on the base class unchanged.

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

* Negate a base form operator through its parent type

Interpolate.__neg__ chose between the expression and the form negation by
testing is_dual(self._function_space), the interpolation's target space. That
space stays primal even when the interpolation is itself a form, so negating an
Interpolate built against a Cofunction gave a Product while multiplying the
same object by -1 gave a FormSum.

None of the other arithmetic operators consult the target space. __add__,
__radd__, __mul__ and __rmul__ all dispatch through _parent_type, which tests
is_dual(self.ufl_function_space()) -- the space the operator takes its value
in, and the one that is dual for a form. Routing __neg__ the same way makes
negation and multiplication by -1 agree by construction.

Firedrake's fml is what notices: it builds `form - label(form)` and compares
the result against `-form`.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_014JYLXVECzcJXgt9gieUbon

* Extract build_coefficient_split from FormData

The mixed-coefficient split dict construction was inlined in
FormData.__init__ and duplicated by tsfc/driver.py::compile_interpolate.
Neither needs Form/Integral machinery for it, so it becomes a standalone
helper in apply_coefficient_split, callable by both.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RFoiWhAUbP53byrgoqxx7H

* Fix interpolation form signature collisions

* Revert "Negate a base form operator through its parent type"

This reverts commit f98822a.

"Fix shape and negation" gave Interpolate both a ufl_shape and a __neg__ keyed
off _function_space, the interpolation target. The target is primal even when
the interpolation carries a test function, so negating one stayed an
expression. f98822a rerouted __neg__ through _parent_type, which keys off
ufl_function_space() -- the adjoint's source dual -- and left ufl_shape on the
target. Negation then returned a FormSum whose ufl_shape does not exist, so any
integrand subtracting an interpolated test function raised AttributeError:

    inner(grad(w) - Interpolate(beta, R), grad(v) - Interpolate(theta, R))

which is the MITC reduction operator applied to both the trial and the test
function. Bisecting that expression across this stack puts the first failure
exactly at f98822a.

_parent_type still drives __add__, __radd__, __mul__ and __rmul__, which
test_interpolate_expr requires: a sum of adjoint interpolations has to stay a
FormSum so it can be passed as the second argument to Interpolate. Only
negation belongs on the target.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JYN9RQAEu7DCRYNkYQGSkW

* Scale a base form operator by the space its value is in

Negation and scalar multiplication disagreed for an interpolation of a test
function: __neg__ followed _function_space, the target, while __rmul__ went
through _parent_type and so followed ufl_function_space(), the adjoint's source
dual. -Interpolate(v, R) was a Product and -1*Interpolate(v, R) a FormSum.

f98822a closed that gap by moving __neg__ onto the form side. That direction
loses the expression reading: FormSum has no ufl_shape, so an integrand that
subtracts an interpolated test function stopped building. Close it from the
other side instead, and let scaling follow the target, which is the space
ufl_shape already reports.

_parent_type still drives __add__ and __radd__, where a sum of adjoint
interpolations has to stay a FormSum for test_interpolate_expr.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JYN9RQAEu7DCRYNkYQGSkW

* Test that scaling an interpolation agrees with negating it

The two disagreed for an interpolation of a test function, and nothing covered
it: test_shape_and_negation only checked negation. Assert that -I, -1*I and
I*-1 give the same type for a trial and a test operand, that scaling keeps the
value shape a FormSum would not have, and that addition still returns a
FormSum, which is the asymmetry test_interpolate_expr depends on.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JYN9RQAEu7DCRYNkYQGSkW

---------

Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
@pbrubeck pbrubeck closed this Sep 11, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants