Fix Interpolate holes

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,)))

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

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,)))

test/test_degree_estimation.py

@@ -24,6 +24,7 @@
     triangle,
 )
 from ufl.algorithms import estimate_total_polynomial_degree
+from ufl.core.interpolate import Interpolate
 
 
 def test_total_degree_estimation():
@@ -58,6 +59,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
 
     # TODO: This should be 1, but 2 is expected behaviour now

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,
@@ -32,8 +34,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 +90,76 @@ 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_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)
+
+    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):

ufl/algorithms/apply_coefficient_split.py

@@ -15,6 +15,7 @@
     Coefficient,
     ComponentTensor,
     Expr,
+    Interpolate,
     MultiIndex,
     NegativeRestricted,
     PositiveRestricted,
@@ -154,6 +155,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,

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)
 

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.
 

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):

ufl/algorithms/estimate_degrees.py

@@ -92,6 +92,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.