From fbdb7d7d88b686c92b38f16ebeaa88835676ed88 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 8 Sep 2026 01:47:32 +0100 Subject: [PATCH 1/2] 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. From 078a1cdc81a5531959e651efaf6e7d03b409ac63 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 11 Sep 2026 08:13:12 +0000 Subject: [PATCH 2/2] Fix Interpolate holes (#497) * 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 * 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 * 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 * 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 * Fix interpolation form signature collisions * 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 * 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 * 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 --------- Co-authored-by: Claude Opus 5 --- 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 | 122 +++++++++++++++++++++- ufl/algorithms/apply_coefficient_split.py | 47 +++++++++ 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/algorithms/formdata.py | 18 ++-- ufl/algorithms/signature.py | 7 ++ ufl/core/base_form_operator.py | 24 ++++- ufl/core/interpolate.py | 119 +++++++++++++++++++++ ufl/differentiation.py | 1 + ufl/domain.py | 10 ++ ufl/referencevalue.py | 9 +- 17 files changed, 395 insertions(+), 23 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..78513718d 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, @@ -14,8 +14,11 @@ Cofunction, FunctionSpace, Mesh, + SpatialCoordinate, TestFunction, + TestFunctions, TrialFunction, + TrialFunctions, action, adjoint, derivative, @@ -32,8 +35,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 Product, ReferenceGrad, ReferenceValue from ufl.core.interpolate import Interpolate from ufl.form import Form, FormSum from ufl.pullback import identity_pullback @@ -85,6 +91,120 @@ 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() + 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_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_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) + 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_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) + + 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..58b43efb3 100644 --- a/ufl/algorithms/apply_coefficient_split.py +++ b/ufl/algorithms/apply_coefficient_split.py @@ -15,6 +15,9 @@ Coefficient, ComponentTensor, Expr, + FormArgument, + FunctionSpace, + Interpolate, MultiIndex, NegativeRestricted, PositiveRestricted, @@ -26,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 @@ -155,6 +159,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, @@ -215,6 +238,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) @@ -266,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/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/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 = [] 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/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 fada73491..4c56901e7 100644 --- a/ufl/core/interpolate.py +++ b/ufl/core/interpolate.py @@ -8,11 +8,17 @@ # # Modified by Nacime Bouziani, 2021-2022 +import hashlib +from collections import defaultdict + 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.operator import Operator 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 @@ -67,12 +73,125 @@ 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 BaseFormOperator.__init__( self, operand, function_space=function_space, argument_slots=argument_slots ) + self._cache = {} + 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 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 = ( + 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 + + 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 _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.""" 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. diff --git a/ufl/domain.py b/ufl/domain.py index 4a3a9328d..750885e60 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] 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