From 1107a482e072bbffa99fa05ce6aa412ac841a83d Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 16 Jul 2026 09:45:28 +0100 Subject: [PATCH 01/16] Fix Interpolate holes --- test/test_apply_coefficent_split.py | 11 +++++++++++ test/test_apply_restrictions.py | 8 ++++++++ test/test_check_arities.py | 12 +++++++++++- test/test_degree_estimation.py | 2 ++ test/test_interpolate.py | 13 +++++++++++++ ufl/algorithms/apply_coefficient_split.py | 20 ++++++++++++++++++++ ufl/algorithms/apply_derivatives.py | 4 ++-- ufl/algorithms/apply_restrictions.py | 11 +++++++++-- ufl/algorithms/check_arities.py | 1 + ufl/algorithms/estimate_degrees.py | 12 ++++++++++++ ufl/core/interpolate.py | 5 +++++ ufl/referencevalue.py | 9 +++++---- 12 files changed, 99 insertions(+), 9 deletions(-) diff --git a/test/test_apply_coefficent_split.py b/test/test_apply_coefficent_split.py index 197bfe570..6120f62ea 100644 --- a/test/test_apply_coefficent_split.py +++ b/test/test_apply_coefficent_split.py @@ -12,6 +12,7 @@ from ufl.classes import ( ComponentTensor, Indexed, + Interpolate, ListTensor, PositiveRestricted, ReferenceGrad, @@ -67,6 +68,16 @@ def test_apply_coefficient_split(self): assert idx1_ == idx1 +def test_interpolate_is_terminal_modifier_boundary(): + cell = triangle + mesh = Mesh(LagrangeElement(cell, 1, (2,))) + V = FunctionSpace(mesh, LagrangeElement(cell, 2)) + interpolation = Interpolate(Coefficient(V), V) + expr = PositiveRestricted(ReferenceGrad(ReferenceValue(interpolation))) + + assert apply_coefficient_split(expr, {}) == expr + + def test_derivative_zero_simplication(): cell = triangle mesh = Mesh(LagrangeElement(cell, 1, (2,))) diff --git a/test/test_apply_restrictions.py b/test/test_apply_restrictions.py index 8a3c84026..70fdd9e4e 100755 --- a/test/test_apply_restrictions.py +++ b/test/test_apply_restrictions.py @@ -14,6 +14,7 @@ ) from ufl.algorithms.apply_restrictions import apply_restrictions from ufl.algorithms.renumbering import renumber_indices +from ufl.core.interpolate import Interpolate from ufl.pullback import identity_pullback from ufl.sobolevspace import L2 @@ -46,6 +47,13 @@ def test_apply_restrictions(): # provided otherwise the user choice is respected assert apply_restrictions(f, default_restrictions={domain: "+"}) == f("+") assert apply_restrictions(f("-"), default_restrictions={domain: "+"}) == f("-") + interpolation = Interpolate(f, v2_space) + assert apply_restrictions(interpolation, default_restrictions={domain: "+"}) == interpolation( + "+" + ) + discontinuous_interpolation = Interpolate(f, v0_space) + with pytest.raises(BaseException): + apply_restrictions(discontinuous_interpolation, default_restrictions={domain: "+"}) assert apply_restrictions(f("+"), default_restrictions={domain: "+"}) == f("+") # Propagation to terminals diff --git a/test/test_check_arities.py b/test/test_check_arities.py index d35eacc9d..b9932b802 100755 --- a/test/test_check_arities.py +++ b/test/test_check_arities.py @@ -21,8 +21,9 @@ inner, tetrahedron, ) -from ufl.algorithms.check_arities import ArityMismatch +from ufl.algorithms.check_arities import ArityMismatch, check_integrand_arity from ufl.algorithms.compute_form_data import compute_form_data +from ufl.core.interpolate import Interpolate def test_check_arities(): @@ -50,6 +51,15 @@ def test_check_arities(): compute_form_data(a) +def test_interpolate_arity(): + domain = Mesh(LagrangeElement(tetrahedron, 1, (3,))) + V = FunctionSpace(domain, LagrangeElement(tetrahedron, 2)) + v = TestFunction(V) + u = TrialFunction(V) + + check_integrand_arity(inner(Interpolate(u, V), v), (v, u)) + + def test_complex_arities(): cell = tetrahedron D = Mesh(LagrangeElement(cell, 1, (3,))) diff --git a/test/test_degree_estimation.py b/test/test_degree_estimation.py index f5fae04c8..93a08edf8 100755 --- a/test/test_degree_estimation.py +++ b/test/test_degree_estimation.py @@ -27,6 +27,7 @@ triangle, ) from ufl.algorithms import estimate_total_polynomial_degree +from ufl.core.interpolate import Interpolate def test_total_degree_estimation(): @@ -61,6 +62,7 @@ def test_total_degree_estimation(): assert estimate_total_polynomial_degree(vu[i] * vv[i]) == 6 assert estimate_total_polynomial_degree(v1) == 1 + assert estimate_total_polynomial_degree(Interpolate(Coefficient(v1_space), v2_space)) == 2 assert estimate_total_polynomial_degree(v2) == 2 # f1 lives on the mixed element's degree-1 sub-element, so its diff --git a/test/test_interpolate.py b/test/test_interpolate.py index d741dc76a..71a81405a 100644 --- a/test/test_interpolate.py +++ b/test/test_interpolate.py @@ -32,8 +32,11 @@ extract_base_form_operators, extract_coefficients, extract_terminals_with_domain, + extract_type, ) +from ufl.algorithms.apply_derivatives import apply_derivatives from ufl.algorithms.expand_indices import expand_indices +from ufl.classes import ReferenceGrad, ReferenceValue from ufl.core.interpolate import Interpolate from ufl.form import Form, FormSum from ufl.pullback import identity_pullback @@ -85,6 +88,16 @@ def test_symbolic(V1, V2): assert Iu.argument_slots() == (vstar, u) assert Iu.arguments() == (vstar,) assert Iu.ufl_operands == (u,) + assert Iu.ufl_element() == V2.ufl_element() + + +def test_reference_value_derivative(V1, V2): + Iu = Interpolate(Coefficient(V1), V2) + reference_value = ReferenceValue(Iu) + + expression = apply_derivatives(grad(reference_value)) + reference_grads = extract_type(expression, ReferenceGrad) + assert reference_value in {g.ufl_operands[0] for g in reference_grads} def test_symbolic_adjoint(V1, V2): diff --git a/ufl/algorithms/apply_coefficient_split.py b/ufl/algorithms/apply_coefficient_split.py index 0bc59cd43..544233700 100644 --- a/ufl/algorithms/apply_coefficient_split.py +++ b/ufl/algorithms/apply_coefficient_split.py @@ -15,6 +15,7 @@ Coefficient, ComponentTensor, Expr, + Interpolate, MultiIndex, NegativeRestricted, PositiveRestricted, @@ -155,6 +156,25 @@ def _( restricted=o._side, ) + @process.register(Interpolate) + def _( + self, + o: Interpolate, + reference_value: bool | None = False, + reference_grad: int = 0, + restricted: str | None = None, + ) -> Expr: + """Handle Interpolate as a finite element terminal.""" + dual_arg, operand = o.argument_slots() + operand = self(operand) + o = o._ufl_expr_reconstruct_(operand, v=dual_arg) + return self._handle_terminal( + o, + reference_value=reference_value, + reference_grad=reference_grad, + restricted=restricted, + ) + @process.register(Terminal) def _( self, diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index e1de2ee06..fe79d0991 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -835,7 +835,7 @@ def _(self, o: ReferenceValue) -> Expr: """Differentiate a reference_value.""" # grad(o) == grad(rv(f)) -> K_ji*rgrad(rv(f))_rj f = o.ufl_operands[0] - if not f._ufl_is_terminal_: + if not (f._ufl_is_terminal_ or isinstance(f, Interpolate)): raise ValueError("ReferenceValue can only wrap a terminal") domain = extract_unique_domain(f, expand_mesh_sequence=False) if isinstance(domain, MeshSequence): @@ -1068,7 +1068,7 @@ def _(self, o: Expr) -> Expr: @process.register(ReferenceValue) def _(self, o: Expr) -> Expr: """Differentiate a reference_value.""" - if not o.ufl_operands[0]._ufl_is_terminal_: + if not (o.ufl_operands[0]._ufl_is_terminal_ or isinstance(o.ufl_operands[0], Interpolate)): raise ValueError("ReferenceValue can only wrap a terminal") return ReferenceGrad(o) diff --git a/ufl/algorithms/apply_restrictions.py b/ufl/algorithms/apply_restrictions.py index d6bb61017..bf7a1bb9b 100644 --- a/ufl/algorithms/apply_restrictions.py +++ b/ufl/algorithms/apply_restrictions.py @@ -15,7 +15,7 @@ from typing import Literal from ufl.algorithms.map_integrands import map_integrand_dags -from ufl.classes import Expr, Restricted +from ufl.classes import Expr, Interpolate, Restricted from ufl.corealg.map_dag import map_expr_dag from ufl.corealg.multifunction import MultiFunction from ufl.domain import Mesh, extract_unique_domain @@ -208,7 +208,7 @@ def variable(self, o, op, label): def reference_value(self, o): """Reference value of something follows same restriction rule as the underlying object.""" (f,) = o.ufl_operands - assert f._ufl_is_terminal_ + assert f._ufl_is_terminal_ or isinstance(f, Interpolate) g = self(f) if isinstance(g, Restricted): side = g.side() @@ -257,6 +257,13 @@ def reference_value(self, o): max_facet_edge_length = _default_restricted facet_origin = _default_restricted # FIXME: Is this valid for quads? + def interpolate(self, o): + """Restrict an interpolated finite element field.""" + if o.ufl_element() in H1: + return self._default_restricted(o) + else: + return self._require_restriction(o) + def coefficient(self, o): """Restrict a coefficient. diff --git a/ufl/algorithms/check_arities.py b/ufl/algorithms/check_arities.py index 2d51bbe1d..b91582ca7 100644 --- a/ufl/algorithms/check_arities.py +++ b/ufl/algorithms/check_arities.py @@ -122,6 +122,7 @@ def linear_operator(self, o, a): grad = linear_operator reference_grad = linear_operator reference_value = linear_operator + interpolate = linear_operator # Conj, is a sesquilinear operator def conj(self, o, a): diff --git a/ufl/algorithms/estimate_degrees.py b/ufl/algorithms/estimate_degrees.py index db6e57d74..37e7c7fa4 100644 --- a/ufl/algorithms/estimate_degrees.py +++ b/ufl/algorithms/estimate_degrees.py @@ -98,6 +98,18 @@ def coefficient(self, v): d = self.default_degree return d + def interpolate(self, v, *ops): + """Apply to interpolate. + + An interpolated field has the polynomial degree of its target element. + """ + e = v.ufl_element() + e = self.element_replace_map.get(e, e) + d = e.embedded_superdegree + if d is None: + d = self.default_degree + return d + def _reduce_degree(self, v, f): """Reduce the estimated degree by one. diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index fada73491..e4e318c79 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -13,6 +13,7 @@ from ufl.core.base_form_operator import BaseFormOperator from ufl.core.ufl_type import ufl_type from ufl.duals import is_dual +from ufl.finiteelement import AbstractFiniteElement from ufl.form import BaseForm from ufl.functionspace import AbstractFunctionSpace @@ -74,6 +75,10 @@ def __init__(self, expr, v): self, operand, function_space=function_space, argument_slots=argument_slots ) + def ufl_element(self) -> AbstractFiniteElement: + """Return the target finite element.""" + return self.ufl_function_space().ufl_element() + def _ufl_expr_reconstruct_(self, expr, v=None, **add_kwargs): """Return a new object of the same type with new operands.""" v = v or self.argument_slots()[0] diff --git a/ufl/referencevalue.py b/ufl/referencevalue.py index 5dd82f4c3..e95ce9ac6 100644 --- a/ufl/referencevalue.py +++ b/ufl/referencevalue.py @@ -5,6 +5,7 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later +from ufl.core.interpolate import Interpolate from ufl.core.operator import Operator from ufl.core.terminal import FormArgument from ufl.core.ufl_type import ufl_type @@ -12,14 +13,14 @@ @ufl_type(num_ops=1, is_index_free=True, is_terminal_modifier=True, is_in_reference_frame=True) class ReferenceValue(Operator): - """Representation of the reference cell value of a form argument.""" + """Representation of the reference cell value of a finite element field.""" __slots__ = () - def __init__(self, f): + def __init__(self, f: FormArgument | Interpolate) -> None: """Initialise.""" - if not isinstance(f, FormArgument): - raise ValueError("Can only take reference value of form arguments.") + if not isinstance(f, FormArgument | Interpolate): + raise ValueError("Can only take reference value of finite element fields.") Operator.__init__(self, (f,)) @property From fe36a1160ae765d73dbe2060c725ff8f209a5551 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 16 Jul 2026 21:35:16 +0100 Subject: [PATCH 02/16] Interpolate._cache --- test/test_interpolate.py | 1 + ufl/core/interpolate.py | 1 + 2 files changed, 2 insertions(+) diff --git a/test/test_interpolate.py b/test/test_interpolate.py index 71a81405a..3ccc5a98b 100644 --- a/test/test_interpolate.py +++ b/test/test_interpolate.py @@ -89,6 +89,7 @@ def test_symbolic(V1, V2): assert Iu.arguments() == (vstar,) assert Iu.ufl_operands == (u,) assert Iu.ufl_element() == V2.ufl_element() + assert Iu._cache == {} def test_reference_value_derivative(V1, V2): diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index e4e318c79..011f9d1e1 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -74,6 +74,7 @@ def __init__(self, expr, v): BaseFormOperator.__init__( self, operand, function_space=function_space, argument_slots=argument_slots ) + self._cache = {} def ufl_element(self) -> AbstractFiniteElement: """Return the target finite element.""" From f8cacc9322102652cb54f9e1e0e6e4461b98058f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 17 Jul 2026 10:25:52 +0100 Subject: [PATCH 03/16] Expose Interpolate to the form compiler --- test/test_interpolate.py | 38 ++++++++++++++ ufl/core/interpolate.py | 110 ++++++++++++++++++++++++++++++++++++++- 2 files changed, 147 insertions(+), 1 deletion(-) diff --git a/test/test_interpolate.py b/test/test_interpolate.py index 3ccc5a98b..9df8b34ca 100644 --- a/test/test_interpolate.py +++ b/test/test_interpolate.py @@ -92,6 +92,44 @@ def test_symbolic(V1, V2): assert Iu._cache == {} +def test_form_compiler_metadata(domain_2d, V1, V2): + u = Coefficient(V1) + cofunction = Cofunction(V2.dual()) + + interpolation = Interpolate(u, V2) + assert interpolation.coefficients() == (u,) + assert interpolation.ufl_domains() == (domain_2d,) + assert interpolation.subdomain_data() == { + domain_2d: {"cell": [None]} + } + assert interpolation.ufl_element() == V2.ufl_element() + + adjoint_interpolation = Interpolate(TestFunction(V1), cofunction) + assert adjoint_interpolation.coefficients() == (cofunction,) + assert adjoint_interpolation.ufl_function_space() == V1.dual() + assert adjoint_interpolation.ufl_element() == V2.ufl_element() + + scalar_interpolation = Interpolate(u, cofunction) + assert scalar_interpolation.arguments() == () + assert scalar_interpolation.coefficients() == (u, cofunction) + assert scalar_interpolation.ufl_function_space() is None + assert scalar_interpolation.ufl_element() == V2.ufl_element() + + +def test_form_compiler_signature(V1, V2, V3): + interpolation = Interpolate(Coefficient(V1), V2) + equivalent = Interpolate(Coefficient(V1), V2) + assert interpolation.signature() == equivalent.signature() + + nested = Interpolate(interpolation, V3) + different_inner_target = Interpolate(Interpolate(Coefficient(V1), V1), V3) + assert nested.signature() != different_inner_target.signature() + + cofunction_sum = Cofunction(V1.dual()) + Cofunction(V1.dual()) + adjoint_interpolation = Interpolate(TestFunction(V2), cofunction_sum) + assert isinstance(adjoint_interpolation.signature(), str) + + def test_reference_value_derivative(V1, V2): Iu = Interpolate(Coefficient(V1), V2) reference_value = ReferenceValue(Iu) diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index 011f9d1e1..184898b6c 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -8,7 +8,12 @@ # # Modified by Nacime Bouziani, 2021-2022 +import hashlib +from collections import defaultdict +from itertools import chain + from ufl.argument import Argument, Coargument +from ufl.coefficient import Cofunction from ufl.constantvalue import as_ufl from ufl.core.base_form_operator import BaseFormOperator from ufl.core.ufl_type import ufl_type @@ -68,6 +73,7 @@ def __init__(self, expr, v): argument_slots = (v, expr) # Get the primal space (V** = V) function_space = v.arguments()[0].ufl_function_space() + self._function_space = function_space # Set the operand as `expr` for DAG traversal purpose. operand = expr @@ -75,10 +81,112 @@ def __init__(self, expr, v): self, operand, function_space=function_space, argument_slots=argument_slots ) self._cache = {} + self._domains = None + self._signature = None + self._subdomain_data = None + self._terminal_numbering = None + + def _analyze_form_arguments(self) -> None: + """Analyze arguments and coefficients in the interpolation.""" + from ufl.algorithms.analysis import extract_coefficients + + super()._analyze_form_arguments() + self._coefficients = tuple(extract_coefficients(self)) + + def _analyze_domains(self) -> None: + """Analyze domains in the interpolation and its argument slots.""" + from ufl.domain import extract_domains, join_domains, sort_domains + + def extract(expression): + if isinstance(expression, BaseForm): + return expression.ufl_domains() + return extract_domains(expression) + + expressions = (*self.ufl_operands, *self.argument_slots()) + self._domains = sort_domains( + join_domains(chain.from_iterable(extract(e) for e in expressions)) + ) + + def ufl_domains(self): + """Return all domains found in the interpolation.""" + if self._domains is None: + self._analyze_domains() + return self._domains + + def subdomain_data(self): + """Return cell-iteration subdomain data for the target domain.""" + if self._subdomain_data is None: + domain = self._function_space.ufl_domain() + self._subdomain_data = {domain: {"cell": [None]}} + return self._subdomain_data + + def terminal_numbering(self): + """Return a contiguous numbering for counted interpolation objects.""" + from ufl.algorithms.analysis import extract_type + from ufl.utils.counted import Counted + from ufl.utils.sorting import sorted_by_count + + if self._terminal_numbering is None: + exprs_by_type = defaultdict(set) + for counted_expr in extract_type(self, Counted): + exprs_by_type[counted_expr._counted_class].add(counted_expr) + + numbering = { + expression: i for i, expression in enumerate(self.arguments()) + } + numbering.update( + { + expression: i + for i, expression in enumerate(self.coefficients()) + } + ) + for expressions in exprs_by_type.values(): + for i, expression in enumerate(sorted_by_count(expressions)): + numbering.setdefault(expression, i) + self._terminal_numbering = numbering + return self._terminal_numbering + + def signature(self): + """Return a numbering-independent signature for compiler caches.""" + from ufl.algorithms.signature import compute_expression_signature + from ufl.form import Form, FormSum + + if self._signature is None: + renumbering = {domain: i for i, domain in enumerate(self.ufl_domains())} + renumbering.update(self.terminal_numbering()) + + def signature(slot): + if isinstance(slot, Interpolate): + return "Interpolate", slot.signature() + if isinstance(slot, Form): + return "Form", slot.signature() + if isinstance(slot, FormSum): + return "FormSum", tuple( + (signature(component), signature(as_ufl(weight))) + for component, weight in zip( + slot.components(), slot.weights() + ) + ) + if isinstance(slot, Coargument | Cofunction): + kind = type(slot).__name__ + slot = Argument(slot.ufl_function_space().dual(), 0) + renumbering[slot] = 0 + return kind, compute_expression_signature(slot, renumbering) + if isinstance(slot, BaseForm): + return type(slot).__name__, tuple( + signature(operand) for operand in slot.ufl_operands + ) + return compute_expression_signature(slot, renumbering) + + signatures = tuple( + signature(slot) for slot in self.argument_slots() + ) + self._signature = hashlib.sha512(str(signatures).encode("utf-8")).hexdigest() + return self._signature def ufl_element(self) -> AbstractFiniteElement: """Return the target finite element.""" - return self.ufl_function_space().ufl_element() + return self._function_space.ufl_element() def _ufl_expr_reconstruct_(self, expr, v=None, **add_kwargs): """Return a new object of the same type with new operands.""" From 9754abc62a0bc4cd2ee0d6b1af89ba27a012b3f8 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 17 Jul 2026 14:45:10 +0100 Subject: [PATCH 04/16] Fix shape and negation --- test/test_interpolate.py | 27 +++++++++++++++++++++++++-- ufl/core/interpolate.py | 13 +++++++++++++ 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/test/test_interpolate.py b/test/test_interpolate.py index 9df8b34ca..2cc3ae62d 100644 --- a/test/test_interpolate.py +++ b/test/test_interpolate.py @@ -4,7 +4,7 @@ __date__ = "2021-11-19" import pytest -from utils import FiniteElement, LagrangeElement +from utils import FiniteElement, LagrangeElement, MixedElement from ufl import ( Action, @@ -15,7 +15,9 @@ FunctionSpace, Mesh, TestFunction, + TestFunctions, TrialFunction, + TrialFunctions, action, adjoint, derivative, @@ -36,7 +38,7 @@ ) from ufl.algorithms.apply_derivatives import apply_derivatives from ufl.algorithms.expand_indices import expand_indices -from ufl.classes import ReferenceGrad, ReferenceValue +from ufl.classes import Product, ReferenceGrad, ReferenceValue from ufl.core.interpolate import Interpolate from ufl.form import Form, FormSum from ufl.pullback import identity_pullback @@ -116,6 +118,27 @@ def test_form_compiler_metadata(domain_2d, V1, V2): assert scalar_interpolation.ufl_element() == V2.ufl_element() +def test_shape_and_negation(domain_2d, V1, V2): + scalar_element = V1.ufl_element() + vector_element = FiniteElement( + "CG", triangle, 1, (2,), identity_pullback, H1 + ) + mixed_space = FunctionSpace( + domain_2d, MixedElement([scalar_element, vector_element]) + ) + target_space = FunctionSpace(domain_2d, vector_element) + _, trial = TrialFunctions(mixed_space) + _, test = TestFunctions(mixed_space) + + for argument in (trial, test): + interpolation = Interpolate(argument, target_space) + assert interpolation.ufl_shape == target_space.value_shape + assert not isinstance(-interpolation, FormSum) + + assert isinstance(-Interpolate(Coefficient(V1), V2), Product) + assert isinstance(-Interpolate(Coefficient(V1), Cofunction(V2.dual())), Product) + + def test_form_compiler_signature(V1, V2, V3): interpolation = Interpolate(Coefficient(V1), V2) equivalent = Interpolate(Coefficient(V1), V2) diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index 184898b6c..7a3ab5044 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -16,6 +16,7 @@ from ufl.coefficient import Cofunction from ufl.constantvalue import as_ufl from ufl.core.base_form_operator import BaseFormOperator +from ufl.core.operator import Operator from ufl.core.ufl_type import ufl_type from ufl.duals import is_dual from ufl.finiteelement import AbstractFiniteElement @@ -188,6 +189,18 @@ def ufl_element(self) -> AbstractFiniteElement: """Return the target finite element.""" return self._function_space.ufl_element() + @property + def ufl_shape(self): + """Return the value shape in the interpolation target space.""" + return self._function_space.value_shape + + def __neg__(self): + """Negate the interpolation result.""" + function_space = self._function_space + if function_space is None or not is_dual(function_space): + return Operator.__rmul__(self, -1) + return BaseForm.__neg__(self) + def _ufl_expr_reconstruct_(self, expr, v=None, **add_kwargs): """Return a new object of the same type with new operands.""" v = v or self.argument_slots()[0] From a324d058181c2da4f0be5c4f5073ceea28da8e6e Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 27 Jul 2026 16:06:52 +0100 Subject: [PATCH 05/16] ruff --- test/test_interpolate.py | 12 +++--------- ufl/core/interpolate.py | 19 ++++--------------- 2 files changed, 7 insertions(+), 24 deletions(-) diff --git a/test/test_interpolate.py b/test/test_interpolate.py index 2cc3ae62d..1b93d1250 100644 --- a/test/test_interpolate.py +++ b/test/test_interpolate.py @@ -101,9 +101,7 @@ def test_form_compiler_metadata(domain_2d, V1, V2): interpolation = Interpolate(u, V2) assert interpolation.coefficients() == (u,) assert interpolation.ufl_domains() == (domain_2d,) - assert interpolation.subdomain_data() == { - domain_2d: {"cell": [None]} - } + assert interpolation.subdomain_data() == {domain_2d: {"cell": [None]}} assert interpolation.ufl_element() == V2.ufl_element() adjoint_interpolation = Interpolate(TestFunction(V1), cofunction) @@ -120,12 +118,8 @@ def test_form_compiler_metadata(domain_2d, V1, V2): def test_shape_and_negation(domain_2d, V1, V2): scalar_element = V1.ufl_element() - vector_element = FiniteElement( - "CG", triangle, 1, (2,), identity_pullback, H1 - ) - mixed_space = FunctionSpace( - domain_2d, MixedElement([scalar_element, vector_element]) - ) + vector_element = FiniteElement("CG", triangle, 1, (2,), identity_pullback, H1) + mixed_space = FunctionSpace(domain_2d, MixedElement([scalar_element, vector_element])) target_space = FunctionSpace(domain_2d, vector_element) _, trial = TrialFunctions(mixed_space) _, test = TestFunctions(mixed_space) diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index 7a3ab5044..add8a0346 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -132,15 +132,8 @@ def terminal_numbering(self): for counted_expr in extract_type(self, Counted): exprs_by_type[counted_expr._counted_class].add(counted_expr) - numbering = { - expression: i for i, expression in enumerate(self.arguments()) - } - numbering.update( - { - expression: i - for i, expression in enumerate(self.coefficients()) - } - ) + numbering = {expression: i for i, expression in enumerate(self.arguments())} + numbering.update({expression: i for i, expression in enumerate(self.coefficients())}) for expressions in exprs_by_type.values(): for i, expression in enumerate(sorted_by_count(expressions)): numbering.setdefault(expression, i) @@ -164,9 +157,7 @@ def signature(slot): if isinstance(slot, FormSum): return "FormSum", tuple( (signature(component), signature(as_ufl(weight))) - for component, weight in zip( - slot.components(), slot.weights() - ) + for component, weight in zip(slot.components(), slot.weights()) ) if isinstance(slot, Coargument | Cofunction): kind = type(slot).__name__ @@ -179,9 +170,7 @@ def signature(slot): ) return compute_expression_signature(slot, renumbering) - signatures = tuple( - signature(slot) for slot in self.argument_slots() - ) + signatures = tuple(signature(slot) for slot in self.argument_slots()) self._signature = hashlib.sha512(str(signatures).encode("utf-8")).hexdigest() return self._signature From 304b80c207393cd0c2a857c17976b3f779f5ffcb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 27 Jul 2026 16:17:53 +0100 Subject: [PATCH 06/16] Fix mypy --- ufl/algorithms/apply_coefficient_split.py | 2 ++ ufl/core/interpolate.py | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/ufl/algorithms/apply_coefficient_split.py b/ufl/algorithms/apply_coefficient_split.py index 544233700..b59172c3c 100644 --- a/ufl/algorithms/apply_coefficient_split.py +++ b/ufl/algorithms/apply_coefficient_split.py @@ -15,6 +15,7 @@ Coefficient, ComponentTensor, Expr, + FormArgument, Interpolate, MultiIndex, NegativeRestricted, @@ -235,6 +236,7 @@ def _handle_terminal( """Wrap terminal as needed.""" c = o if reference_value: + assert isinstance(c, FormArgument | Interpolate) c = ReferenceValue(c) for k in range(reference_grad): c = ReferenceGrad(c) diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index add8a0346..d53092b9f 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -105,7 +105,7 @@ def extract(expression): expressions = (*self.ufl_operands, *self.argument_slots()) self._domains = sort_domains( - join_domains(chain.from_iterable(extract(e) for e in expressions)) + join_domains(list(chain.from_iterable(extract(e) for e in expressions))) ) def ufl_domains(self): From 4482d080cd709c0351905401c06421ea91d63b64 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 26 Aug 2026 18:12:46 +0100 Subject: [PATCH 07/16] 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 --- ufl/domain.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/ufl/domain.py b/ufl/domain.py index 4a3a9328d..95f07e16d 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -404,9 +404,13 @@ def extract_domains( """ from ufl.algorithms.traversal import iter_expressions + from ufl.core.base_form_operator import BaseFormOperator from ufl.form import Form from ufl.integral import Integral + if isinstance(expr, BaseFormOperator) and expand_mesh_sequence: + # The argument slots carry domains that the operands cannot reach. + return tuple(expr.ufl_domains()) if isinstance(expr, Form): if not expand_mesh_sequence: raise NotImplementedError(""" @@ -450,6 +454,12 @@ def extract_unique_domain( domain. """ + from ufl.core.base_form_operator import BaseFormOperator + + if isinstance(expr, BaseFormOperator): + # A base form operator has the domains of its argument slots as well, + # but it takes its value on the one its operands are defined over. + expr, = expr.ufl_operands domains = extract_domains(expr, expand_mesh_sequence=expand_mesh_sequence) if len(domains) == 1: return domains[0] From ce071590de8a594bbe926e835b89bd54ce8a1e12 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 27 Aug 2026 11:47:07 +0100 Subject: [PATCH 08/16] ruff --- ufl/domain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ufl/domain.py b/ufl/domain.py index 95f07e16d..750885e60 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -459,7 +459,7 @@ def extract_unique_domain( if isinstance(expr, BaseFormOperator): # A base form operator has the domains of its argument slots as well, # but it takes its value on the one its operands are defined over. - expr, = expr.ufl_operands + (expr,) = expr.ufl_operands domains = extract_domains(expr, expand_mesh_sequence=expand_mesh_sequence) if len(domains) == 1: return domains[0] From 9058692d03e6f347e67300b1555bbf28f1af625d Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 27 Aug 2026 12:55:47 +0100 Subject: [PATCH 09/16] 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 --- ufl/core/base_form_operator.py | 24 +++++++++++++++++++++++- ufl/core/interpolate.py | 22 ---------------------- ufl/differentiation.py | 1 + 3 files changed, 24 insertions(+), 23 deletions(-) diff --git a/ufl/core/base_form_operator.py b/ufl/core/base_form_operator.py index edeeb09e5..b7fe6ffb7 100644 --- a/ufl/core/base_form_operator.py +++ b/ufl/core/base_form_operator.py @@ -14,6 +14,7 @@ # Modified by Nacime Bouziani, 2021-2022 from collections import OrderedDict +from itertools import chain from numbers import Number from ufl.argument import Argument, Coargument @@ -70,8 +71,9 @@ def __init__(self, *operands, function_space, derivatives=None, argument_slots=( argument_slots = (v_star,) self._argument_slots = argument_slots - # Internal variables for caching coefficient data + # Internal variables for caching coefficient and domain data self._coefficients = None + self._domains = None # BaseFormOperators don't have free indices. ufl_free_indices = () @@ -96,6 +98,26 @@ def argument_slots(self, outer_form=False): # => F.arguments() should return (v,)! return tuple(a for a in self._argument_slots[1:] if len(extract_arguments(a)) != 0) + def _analyze_domains(self) -> None: + """Analyze domains in the operands and the argument slots.""" + from ufl.domain import extract_domains, join_domains, sort_domains + + def extract(expression): + if isinstance(expression, BaseForm): + return expression.ufl_domains() + return extract_domains(expression) + + expressions = (*self.ufl_operands, *self.argument_slots()) + self._domains = sort_domains( + join_domains(list(chain.from_iterable(extract(e) for e in expressions))) + ) + + def ufl_domains(self): + """Return all domains found in the operands and the argument slots.""" + if self._domains is None: + self._analyze_domains() + return self._domains + def coefficients(self): """Return all BaseCoefficient objects found in base form operator.""" if self._coefficients is None: diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index d53092b9f..2ca827ec6 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -10,7 +10,6 @@ import hashlib from collections import defaultdict -from itertools import chain from ufl.argument import Argument, Coargument from ufl.coefficient import Cofunction @@ -82,7 +81,6 @@ def __init__(self, expr, v): self, operand, function_space=function_space, argument_slots=argument_slots ) self._cache = {} - self._domains = None self._signature = None self._subdomain_data = None self._terminal_numbering = None @@ -94,26 +92,6 @@ def _analyze_form_arguments(self) -> None: super()._analyze_form_arguments() self._coefficients = tuple(extract_coefficients(self)) - def _analyze_domains(self) -> None: - """Analyze domains in the interpolation and its argument slots.""" - from ufl.domain import extract_domains, join_domains, sort_domains - - def extract(expression): - if isinstance(expression, BaseForm): - return expression.ufl_domains() - return extract_domains(expression) - - expressions = (*self.ufl_operands, *self.argument_slots()) - self._domains = sort_domains( - join_domains(list(chain.from_iterable(extract(e) for e in expressions))) - ) - - def ufl_domains(self): - """Return all domains found in the interpolation.""" - if self._domains is None: - self._analyze_domains() - return self._domains - def subdomain_data(self): """Return cell-iteration subdomain data for the target domain.""" if self._subdomain_data is None: diff --git a/ufl/differentiation.py b/ufl/differentiation.py index f660840a3..adf3fbc89 100644 --- a/ufl/differentiation.py +++ b/ufl/differentiation.py @@ -156,6 +156,7 @@ def __init__(self, base_form, coefficients, arguments, coefficient_derivatives): self, base_form, coefficients, arguments, coefficient_derivatives ) self._argument_slots = base_form._argument_slots + self._domains = None # Enforce Operator reconstruction as Operator is a parent class of # both: BaseFormDerivative and BaseFormOperator. From f98822a63696da34a8ec8b14b1e37b5d552c27ab Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 30 Aug 2026 06:41:45 +0100 Subject: [PATCH 10/16] 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 Claude-Session: https://claude.ai/code/session_014JYLXVECzcJXgt9gieUbon --- test/test_interpolate.py | 7 ++++++- ufl/core/interpolate.py | 8 ++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/test/test_interpolate.py b/test/test_interpolate.py index 1b93d1250..ed6f196f6 100644 --- a/test/test_interpolate.py +++ b/test/test_interpolate.py @@ -127,7 +127,12 @@ def test_shape_and_negation(domain_2d, V1, V2): for argument in (trial, test): interpolation = Interpolate(argument, target_space) assert interpolation.ufl_shape == target_space.value_shape - assert not isinstance(-interpolation, FormSum) + # Negation must agree with multiplication by -1, which negates a + # primal interpolation as an expression and a dual one as a form. + assert type(-interpolation) is type(-1 * interpolation) + + assert not isinstance(-Interpolate(trial, target_space), FormSum) + assert isinstance(-Interpolate(test, target_space), FormSum) assert isinstance(-Interpolate(Coefficient(V1), V2), Product) assert isinstance(-Interpolate(Coefficient(V1), Cofunction(V2.dual())), Product) diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index 2ca827ec6..a90a3a348 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -15,7 +15,6 @@ from ufl.coefficient import Cofunction from ufl.constantvalue import as_ufl from ufl.core.base_form_operator import BaseFormOperator -from ufl.core.operator import Operator from ufl.core.ufl_type import ufl_type from ufl.duals import is_dual from ufl.finiteelement import AbstractFiniteElement @@ -162,11 +161,8 @@ def ufl_shape(self): return self._function_space.value_shape def __neg__(self): - """Negate the interpolation result.""" - function_space = self._function_space - if function_space is None or not is_dual(function_space): - return Operator.__rmul__(self, -1) - return BaseForm.__neg__(self) + """Negate.""" + return self._parent_type.__neg__(self) def _ufl_expr_reconstruct_(self, expr, v=None, **add_kwargs): """Return a new object of the same type with new operands.""" From a75e36a1e040b694f8efa91e140581f397fe7cd6 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 7 Sep 2026 00:40:32 +0100 Subject: [PATCH 11/16] 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 Claude-Session: https://claude.ai/code/session_01RFoiWhAUbP53byrgoqxx7H --- ufl/algorithms/apply_coefficient_split.py | 25 +++++++++++++++++++++++ ufl/algorithms/formdata.py | 18 ++++++---------- 2 files changed, 31 insertions(+), 12 deletions(-) diff --git a/ufl/algorithms/apply_coefficient_split.py b/ufl/algorithms/apply_coefficient_split.py index b59172c3c..58b43efb3 100644 --- a/ufl/algorithms/apply_coefficient_split.py +++ b/ufl/algorithms/apply_coefficient_split.py @@ -16,6 +16,7 @@ ComponentTensor, Expr, FormArgument, + FunctionSpace, Interpolate, MultiIndex, NegativeRestricted, @@ -28,6 +29,7 @@ ) from ufl.core.multiindex import indices from ufl.corealg.dag_traverser import DAGTraverser +from ufl.domain import extract_unique_domain from ufl.form import BaseForm from ufl.tensors import as_tensor @@ -288,3 +290,26 @@ def apply_coefficient_split(expr: Expr, coefficient_split: dict) -> Expr: if not coefficient_split: return expr return CoefficientSplitter(coefficient_split)(expr) + + +def build_coefficient_split(coefficients_to_split) -> dict: + """Map each mixed coefficient in ``coefficients_to_split`` to its per-mesh components. + + Args: + coefficients_to_split: Coefficients with a mixed element to split. + + Returns: + `dict` that maps each coefficient to its components, suitable for + `CoefficientSplitter`/`apply_coefficient_split`. + + """ + coefficient_split = {} + for c in coefficients_to_split: + mesh = extract_unique_domain(c, expand_mesh_sequence=False) + assert mesh is not None + elem = c.ufl_element() + coefficient_split[c] = [ + Coefficient(FunctionSpace(m, e)) + for m, e in zip(mesh.iterable_like(elem), elem.sub_elements) + ] + return coefficient_split diff --git a/ufl/algorithms/formdata.py b/ufl/algorithms/formdata.py index 2046b103e..967585606 100644 --- a/ufl/algorithms/formdata.py +++ b/ufl/algorithms/formdata.py @@ -14,7 +14,7 @@ from typing import Any from ufl.algorithms.analysis import extract_coefficients, extract_sub_elements, unique_tuple -from ufl.algorithms.apply_coefficient_split import CoefficientSplitter +from ufl.algorithms.apply_coefficient_split import CoefficientSplitter, build_coefficient_split from ufl.algorithms.apply_restrictions import apply_restrictions, default_restriction_map from ufl.algorithms.check_arities import check_integrand_arity from ufl.algorithms.domain_analysis import IntegralData, reconstruct_form_from_integral_data @@ -278,17 +278,11 @@ def __init__( # Split coefficients that are contained in ``coefficients_to_split`` # into components, and store a dict in ``self`` that maps # each coefficient to its components. - coefficient_split = {} - for o in self.reduced_coefficients: - if o in coefficients_to_split: - c = self.function_replace_map[o] - mesh = extract_unique_domain(c, expand_mesh_sequence=False) - elem = c.ufl_element() - coefficient_split[c] = [ - Coefficient(FunctionSpace(m, e)) - for m, e in zip(mesh.iterable_like(elem), elem.sub_elements) # type: ignore - ] - self._coefficient_split = coefficient_split + self._coefficient_split = build_coefficient_split( + self.function_replace_map[o] + for o in self.reduced_coefficients + if o in coefficients_to_split + ) coeff_splitter = CoefficientSplitter(self.coefficient_split) for itg_data in self.integral_data: new_integrals = [] From 7664655df6370d61281746c1bbfba53d478c0070 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 9 Sep 2026 22:40:37 +0100 Subject: [PATCH 12/16] Fix interpolation form signature collisions --- test/test_interpolate.py | 15 +++++++++++++++ ufl/algorithms/signature.py | 7 +++++++ ufl/core/interpolate.py | 5 ++++- 3 files changed, 26 insertions(+), 1 deletion(-) diff --git a/test/test_interpolate.py b/test/test_interpolate.py index ed6f196f6..0ee02e6dc 100644 --- a/test/test_interpolate.py +++ b/test/test_interpolate.py @@ -14,6 +14,7 @@ Cofunction, FunctionSpace, Mesh, + SpatialCoordinate, TestFunction, TestFunctions, TrialFunction, @@ -152,6 +153,20 @@ def test_form_compiler_signature(V1, V2, V3): assert isinstance(adjoint_interpolation.signature(), str) +def test_form_compiler_signature_depends_on_interpolation_target(domain_2d): + """Different target elements must not share an interpolation signature.""" + curl_element = FiniteElement("N1curl", triangle, 1, (2,), identity_pullback, H1) + lagrange_element = FiniteElement("Lagrange", triangle, 1, (2,), identity_pullback, H1) + curl_space = FunctionSpace(domain_2d, curl_element) + lagrange_space = FunctionSpace(domain_2d, lagrange_element) + x = SpatialCoordinate(domain_2d) + + curl_form = Interpolate(x, curl_space)[0] * dx + lagrange_form = Interpolate(x, lagrange_space)[0] * dx + + assert curl_form.signature() != lagrange_form.signature() + + def test_reference_value_derivative(V1, V2): Iu = Interpolate(Coefficient(V1), V2) reference_value = ReferenceValue(Iu) diff --git a/ufl/algorithms/signature.py b/ufl/algorithms/signature.py index ced238f9c..f2658b7c6 100644 --- a/ufl/algorithms/signature.py +++ b/ufl/algorithms/signature.py @@ -106,6 +106,13 @@ def compute_expression_hashdata(expression, terminal_hashdata) -> bytes: else: data = [expr._ufl_typecode_] + # Interpolate's target space is an argument slot rather than a UFL + # operand and affects the form signature. + from ufl.core.interpolate import Interpolate + + if isinstance(expr, Interpolate): + data.append(expr.signature()) + for op in expr.ufl_operands: data += [cache[op]] cache[expr] = hashlib.sha512(str(data).encode("utf-8")).digest() diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index a90a3a348..38dc07b1f 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -147,7 +147,10 @@ def signature(slot): ) return compute_expression_signature(slot, renumbering) - signatures = tuple(signature(slot) for slot in self.argument_slots()) + signatures = ( + repr(self.ufl_element()), + *(signature(slot) for slot in self.argument_slots()), + ) self._signature = hashlib.sha512(str(signatures).encode("utf-8")).hexdigest() return self._signature From d9fa05906fff16d59cbb26a0cb9ab6f233102c58 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 11 Sep 2026 00:31:26 +0100 Subject: [PATCH 13/16] Revert "Negate a base form operator through its parent type" This reverts commit f98822a63696da34a8ec8b14b1e37b5d552c27ab. "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. f98822a6 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 f98822a6. _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 Claude-Session: https://claude.ai/code/session_01JYN9RQAEu7DCRYNkYQGSkW --- test/test_interpolate.py | 7 +------ ufl/core/interpolate.py | 8 ++++++-- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/test/test_interpolate.py b/test/test_interpolate.py index 0ee02e6dc..ced2137f3 100644 --- a/test/test_interpolate.py +++ b/test/test_interpolate.py @@ -128,12 +128,7 @@ def test_shape_and_negation(domain_2d, V1, V2): for argument in (trial, test): interpolation = Interpolate(argument, target_space) assert interpolation.ufl_shape == target_space.value_shape - # Negation must agree with multiplication by -1, which negates a - # primal interpolation as an expression and a dual one as a form. - assert type(-interpolation) is type(-1 * interpolation) - - assert not isinstance(-Interpolate(trial, target_space), FormSum) - assert isinstance(-Interpolate(test, target_space), FormSum) + assert not isinstance(-interpolation, FormSum) assert isinstance(-Interpolate(Coefficient(V1), V2), Product) assert isinstance(-Interpolate(Coefficient(V1), Cofunction(V2.dual())), Product) diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index 38dc07b1f..bd7ea9289 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -15,6 +15,7 @@ from ufl.coefficient import Cofunction from ufl.constantvalue import as_ufl from ufl.core.base_form_operator import BaseFormOperator +from ufl.core.operator import Operator from ufl.core.ufl_type import ufl_type from ufl.duals import is_dual from ufl.finiteelement import AbstractFiniteElement @@ -164,8 +165,11 @@ def ufl_shape(self): return self._function_space.value_shape def __neg__(self): - """Negate.""" - return self._parent_type.__neg__(self) + """Negate the interpolation result.""" + function_space = self._function_space + if function_space is None or not is_dual(function_space): + return Operator.__rmul__(self, -1) + return BaseForm.__neg__(self) def _ufl_expr_reconstruct_(self, expr, v=None, **add_kwargs): """Return a new object of the same type with new operands.""" From 7ea8b00d76cd46a7f25a405474666650f6ca69aa Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 11 Sep 2026 00:35:08 +0100 Subject: [PATCH 14/16] 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. f98822a6 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 Claude-Session: https://claude.ai/code/session_01JYN9RQAEu7DCRYNkYQGSkW --- ufl/core/interpolate.py | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py index bd7ea9289..4c56901e7 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -164,13 +164,35 @@ def ufl_shape(self): """Return the value shape in the interpolation target space.""" return self._function_space.value_shape - def __neg__(self): - """Negate the interpolation result.""" + def _value_parent_type(self): + """Return the type whose arithmetic matches the interpolation's value. + + An interpolation takes its value in the target space, which is the + space ``ufl_shape`` reports. Scaling and negation follow that space, so + that an interpolation of a test function stays an expression and can be + combined with one inside an integrand. ``_parent_type`` instead follows + ``ufl_function_space()``, the adjoint's source dual, and still drives + addition, where a sum of adjoint interpolations must stay a form. + """ function_space = self._function_space if function_space is None or not is_dual(function_space): + return Operator + return BaseForm + + def __neg__(self): + """Negate the interpolation result.""" + if self._value_parent_type() is Operator: return Operator.__rmul__(self, -1) return BaseForm.__neg__(self) + def __mul__(self, other): + """Multiply, agreeing with negation on which space the value is in.""" + return self._value_parent_type().__mul__(self, other) + + def __rmul__(self, other): + """Multiply, agreeing with negation on which space the value is in.""" + return self._value_parent_type().__rmul__(self, other) + def _ufl_expr_reconstruct_(self, expr, v=None, **add_kwargs): """Return a new object of the same type with new operands.""" v = v or self.argument_slots()[0] From b3cb894daee5974f595455e642bdd94a239c4588 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 11 Sep 2026 00:37:41 +0100 Subject: [PATCH 15/16] 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 Claude-Session: https://claude.ai/code/session_01JYN9RQAEu7DCRYNkYQGSkW --- test/test_interpolate.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/test/test_interpolate.py b/test/test_interpolate.py index ced2137f3..78513718d 100644 --- a/test/test_interpolate.py +++ b/test/test_interpolate.py @@ -134,6 +134,42 @@ def test_shape_and_negation(domain_2d, V1, V2): assert isinstance(-Interpolate(Coefficient(V1), Cofunction(V2.dual())), Product) +def test_scaling_agrees_with_negation(domain_2d, V1, V2): + """Scaling an interpolation must agree with negating it. + + Both follow the target space, the one that ``ufl_shape`` reports, so that an + interpolation of a test function stays an expression. Firedrake's fml is + what notices a disagreement: it builds ``form - label(form)`` and compares + the result against ``-form``. + """ + scalar_element = V1.ufl_element() + vector_element = FiniteElement("CG", triangle, 1, (2,), identity_pullback, H1) + mixed_space = FunctionSpace(domain_2d, MixedElement([scalar_element, vector_element])) + target_space = FunctionSpace(domain_2d, vector_element) + _, trial = TrialFunctions(mixed_space) + _, test = TestFunctions(mixed_space) + + for argument in (trial, test): + interpolation = Interpolate(argument, target_space) + assert type(-interpolation) is type(-1 * interpolation) + assert type(-interpolation) is type(interpolation * -1) + assert not isinstance(2 * interpolation, FormSum) + # Scaling preserves the value shape, which a FormSum would not have. + assert (2 * interpolation).ufl_shape == target_space.value_shape + + for interpolation in ( + Interpolate(Coefficient(V1), V2), + Interpolate(Coefficient(V1), Cofunction(V2.dual())), + ): + assert isinstance(-1 * interpolation, Product) + assert type(-interpolation) is type(-1 * interpolation) + + # Addition still follows ufl_function_space(), so that a sum of adjoint + # interpolations stays a form -- see test_interpolate_expr. + adjoint_interpolation = Interpolate(test, target_space) + assert isinstance(adjoint_interpolation + adjoint_interpolation, FormSum) + + def test_form_compiler_signature(V1, V2, V3): interpolation = Interpolate(Coefficient(V1), V2) equivalent = Interpolate(Coefficient(V1), V2) From fbdb7d7d88b686c92b38f16ebeaa88835676ed88 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 8 Sep 2026 01:47:32 +0100 Subject: [PATCH 16/16] Give each pull back an inverse 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 Claude-Session: https://claude.ai/code/session_0155wQMTMSGiF9a6TbyTVK5P --- test/test_pullback_inverse.py | 98 ++++++++++++++++ ufl/pullback.py | 207 ++++++++++++++++++++++++++++++++++ 2 files changed, 305 insertions(+) create mode 100644 test/test_pullback_inverse.py diff --git a/test/test_pullback_inverse.py b/test/test_pullback_inverse.py new file mode 100644 index 000000000..80c7cd2dd --- /dev/null +++ b/test/test_pullback_inverse.py @@ -0,0 +1,98 @@ +"""Tests of the inverse pull backs.""" + +import numpy as np +import pytest +from utils import FiniteElement, LagrangeElement, MixedElement, SymmetricElement + +from ufl import Cell, Coefficient, FunctionSpace, Mesh +from ufl.algorithms.apply_derivatives import apply_derivatives +from ufl.algorithms.cancel_jacobian_products import cancel_jacobian_products +from ufl.algorithms.remove_component_tensors import remove_component_tensors +from ufl.algorithms.renumbering import renumber_indices +from ufl.classes import JacobianDeterminant, ReferenceValue +from ufl.pullback import ( + contravariant_piola, + covariant_contravariant_piola, + covariant_piola, + double_contravariant_piola, + double_covariant_piola, + l2_piola, + physical_pullback, + undefined_pullback, +) +from ufl.sobolevspace import H1, L2, HCurl, HDiv, HDivDiv, HEin + +cell = Cell("triangle") +domain = Mesh(LagrangeElement(cell, 1, (2,))) + +U = LagrangeElement(cell, 1) +Vd = FiniteElement("Raviart-Thomas", cell, 1, (2,), contravariant_piola, HDiv) +Vc = FiniteElement("N1curl", cell, 1, (2,), covariant_piola, HCurl) +Td = FiniteElement("Regge", cell, 1, (2, 2), double_covariant_piola, HEin) +Tc = FiniteElement("HHJ", cell, 1, (2, 2), double_contravariant_piola, HDivDiv) +Tcc = FiniteElement("CC", cell, 1, (2, 2), covariant_contravariant_piola, HDivDiv) +S = SymmetricElement({(0, 0): 0, (1, 0): 1, (0, 1): 1, (1, 1): 2}, [U, U, U]) +M = MixedElement([U, Vd, Vc]) + + +def simplify(expr): + """Cancel the Jacobian products a round trip leaves behind.""" + return cancel_jacobian_products(remove_component_tensors(apply_derivatives(expr))) + + +@pytest.mark.parametrize( + "element", + [U, Vd, Vc, Td, Tc, Tcc, S, M], + ids=[ + "identity", + "contravariant", + "covariant", + "double covariant", + "double contravariant", + "covariant contravariant", + "symmetric", + "mixed", + ], +) +def test_apply_inverse_undoes_apply(element): + """The inverse pull back returns a pushed-forward function unchanged.""" + pullback = element.pullback + reference = ReferenceValue(Coefficient(FunctionSpace(domain, element))) + actual = simplify(pullback.apply_inverse(pullback.apply(reference))) + + assert actual.ufl_shape == reference.ufl_shape + for idx in np.ndindex(reference.ufl_shape): + assert renumber_indices(actual[idx]) == renumber_indices(reference[idx]) + + +def test_l2_piola_apply_inverse(): + """The L2 Piola scales by the Jacobian determinant. + + The round trip leaves ``detJ / detJ`` standing, because cancelling a scalar + factor is not something ``cancel_jacobian_products`` does. + """ + element = FiniteElement("Discontinuous Lagrange", cell, 1, (), l2_piola, L2) + reference = ReferenceValue(Coefficient(FunctionSpace(domain, element))) + + assert l2_piola.apply_inverse(reference) == reference * JacobianDeterminant(domain) + + +def test_apply_inverse_of_physical_value_shape(): + """The inverse pull back maps a physical shape to a reference shape.""" + for element in [U, Vd, Vc, Td, Tc, Tcc, S, M]: + pullback = element.pullback + physical = pullback.apply(ReferenceValue(Coefficient(FunctionSpace(domain, element)))) + assert physical.ufl_shape == pullback.physical_value_shape(element, domain) + assert pullback.apply_inverse(physical).ufl_shape == element.reference_value_shape + + +@pytest.mark.parametrize("pullback", [physical_pullback, undefined_pullback]) +def test_apply_inverse_is_not_defined(pullback): + """A pull back with no standard inverse says so.""" + element = FiniteElement("Custom", cell, 1, (), pullback, H1) + reference = ReferenceValue(Coefficient(FunctionSpace(domain, element))) + if pullback is undefined_pullback: + with pytest.raises(BaseException): + pullback.apply_inverse(reference) + else: + assert pullback.apply_inverse(reference) == reference diff --git a/ufl/pullback.py b/ufl/pullback.py index 8fd7111c7..5c959eb6a 100644 --- a/ufl/pullback.py +++ b/ufl/pullback.py @@ -78,6 +78,17 @@ def apply(self, expr: Expr, domain: AbstractDomain | None = None) -> Expr: """ raise NonStandardPullbackException() + def apply_inverse(self, expr: Expr, domain: AbstractDomain | None = None) -> Expr: + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + raise NonStandardPullbackException() + class IdentityPullback(AbstractPullback): """The identity pull back.""" @@ -102,6 +113,17 @@ def apply(self, expr, domain=None): """ return expr + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + return expr + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain. @@ -147,6 +169,25 @@ def apply(self, expr, domain=None): kj = (*k, j) return as_tensor(transform[i, j] * expr[kj], (*k, i)) + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + from ufl.classes import Jacobian, JacobianDeterminant, JacobianInverse + + domain = domain or extract_unique_domain(expr) + J = Jacobian(domain) + detJ = JacobianDeterminant(J) + K = JacobianInverse(domain) + *k, i, j = indices(len(expr.ufl_shape) + 1) + kj = (*k, j) + return as_tensor(detJ * K[i, j] * expr[kj], (*k, i)) + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain. @@ -191,6 +232,23 @@ def apply(self, expr, domain=None): kj = (*k, j) return as_tensor(K[j, i] * expr[kj], (*k, i)) + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + from ufl.classes import Jacobian + + domain = domain or extract_unique_domain(expr) + J = Jacobian(domain) + *k, i, j = indices(len(expr.ufl_shape) + 1) + kj = (*k, j) + return as_tensor(J[j, i] * expr[kj], (*k, i)) + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain. @@ -232,6 +290,21 @@ def apply(self, expr, domain=None): detJ = JacobianDeterminant(domain) return expr / detJ + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + from ufl.classes import JacobianDeterminant + + domain = domain or extract_unique_domain(expr) + detJ = JacobianDeterminant(domain) + return expr * detJ + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain. @@ -276,6 +349,24 @@ def apply(self, expr, domain=None): kmn = (*k, m, n) return as_tensor((1.0 / detJ) ** 2 * J[i, m] * expr[kmn] * J[j, n], (*k, i, j)) + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + from ufl.classes import JacobianDeterminant, JacobianInverse + + domain = domain or extract_unique_domain(expr) + detJ = JacobianDeterminant(domain) + K = JacobianInverse(domain) + *k, i, j, m, n = indices(len(expr.ufl_shape) + 2) + kmn = (*k, m, n) + return as_tensor(detJ**2 * K[i, m] * expr[kmn] * K[j, n], (*k, i, j)) + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain. @@ -320,6 +411,23 @@ def apply(self, expr, domain=None): kmn = (*k, m, n) return as_tensor(K[m, i] * expr[kmn] * K[n, j], (*k, i, j)) + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + from ufl.classes import Jacobian + + domain = domain or extract_unique_domain(expr) + J = Jacobian(domain) + *k, i, j, m, n = indices(len(expr.ufl_shape) + 2) + kmn = (*k, m, n) + return as_tensor(J[m, i] * expr[kmn] * J[n, j], (*k, i, j)) + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain. @@ -366,6 +474,25 @@ def apply(self, expr, domain=None): kmn = (*k, m, n) return as_tensor((1.0 / detJ) * K[m, i] * expr[kmn] * J[j, n], (*k, i, j)) + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + from ufl.classes import Jacobian, JacobianDeterminant, JacobianInverse + + domain = domain or extract_unique_domain(expr) + J = Jacobian(domain) + detJ = JacobianDeterminant(J) + K = JacobianInverse(domain) + *k, i, j, m, n = indices(len(expr.ufl_shape) + 2) + kmn = (*k, m, n) + return as_tensor(detJ * J[m, i] * expr[kmn] * K[j, n], (*k, i, j)) + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element. @@ -438,6 +565,33 @@ def apply(self, expr, domain=None): ) return f + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + domain = domain or extract_unique_domain(expr, expand_mesh_sequence=False) + gflat = [expr[idx] for idx in np.ndindex(expr.ufl_shape)] + r_components = [] + offset = 0 + # For each piece in physical space, apply the appropriate inverse pullback + for subelem, subdomain in zip( + self._element.sub_elements, domain.iterable_like(self._element) + ): + physical_shape = subelem.pullback.physical_value_shape(subelem, subdomain) + size = int(np.prod(physical_shape, dtype=int)) + gsub = as_tensor(np.asarray(gflat[offset : offset + size]).reshape(physical_shape)) + gmapped = subelem.pullback.apply_inverse(gsub, domain=subdomain) + # Flatten into the mapped expression for the whole thing + r_components.extend(gmapped[idx] for idx in np.ndindex(gmapped.ufl_shape)) + offset += size + # And reshape appropriately + return as_tensor(np.asarray(r_components).reshape(self._element.reference_value_shape)) + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain. @@ -527,6 +681,37 @@ def apply(self, expr, domain=None): ) return f + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + domain = domain or extract_unique_domain(expr, expand_mesh_sequence=False) + subelem = self._element.sub_elements[0] + physical_shape = subelem.pullback.physical_value_shape(subelem, domain) + size = int(np.prod(physical_shape, dtype=int)) + gflat = [expr[idx] for idx in np.ndindex(expr.ufl_shape)] + # Symmetry repeats a reference piece across several physical blocks, so + # map each piece once, from the first block that carries it. + r_pieces = {} + for block, component in enumerate(np.ndindex(self._block_shape)): + i = self._symmetry[component] + if i in r_pieces: + continue + gsub = as_tensor( + np.asarray(gflat[size * block : size * (block + 1)]).reshape(physical_shape) + ) + r_pieces[i] = subelem.pullback.apply_inverse(gsub, domain=domain) + r_components = [] + for i in sorted(r_pieces): + gmapped = r_pieces[i] + r_components.extend(gmapped[idx] for idx in np.ndindex(gmapped.ufl_shape)) + return as_tensor(np.asarray(r_components).reshape(self._element.reference_value_shape)) + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain. @@ -569,6 +754,17 @@ def apply(self, expr, domain=None): """ return expr + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + return expr + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain. @@ -608,6 +804,17 @@ def apply(self, expr, domain=None): """ return expr + def apply_inverse(self, expr, domain=None): + """Apply the inverse of the pull back. + + Args: + expr: A function on a physical cell + domain: The domain on which the function is defined + + Returns: The function mapped to the reference cell + """ + return expr + def physical_value_shape(self, element, domain) -> tuple[int, ...]: """Get the physical value shape when this pull back is applied to an element on a domain.