interpolate two-forms (#4770)

Author: leo-collins

firedrake/assemble.py

@@ -855,6 +855,15 @@ def restructure_base_form(expr, visited=None):
         if isinstance(expr, ufl.FormSum) and all(ufl.duals.is_dual(a.function_space()) for a in expr.arguments()):
             # Return ufl.Sum if we are assembling a FormSum with Coarguments (a primal expression)
             return sum(w*c for w, c in zip(expr.weights(), expr.components()))
+
+        # If F: V3 x V2 -> R, then
+        # Interpolate(TestFunction(V1), F) <=> Action(Interpolate(TestFunction(V1), TrialFunction(V2.dual())), F).
+        # The result is a two-form V3 x V1 -> R.
+        if isinstance(expr, ufl.Interpolate) and isinstance(expr.argument_slots()[0], ufl.form.Form):
+            form, operand = expr.argument_slots()
+            vstar = firedrake.Argument(form.arguments()[0].function_space().dual(), 1)
+            expr = expr._ufl_expr_reconstruct_(operand, v=vstar)
+            return ufl.action(expr, form)
         return expr
 
     @staticmethod

tests/firedrake/regression/test_interpolate.py

@@ -6,6 +6,13 @@
 cwd = abspath(dirname(__file__))
 
 
+def mat_equals(a, b):
+    """Check that two Matrices are equal."""
+    a = a.petscmat.copy()
+    a.axpy(-1.0, b.petscmat)
+    return a.norm(norm_type=PETSc.NormType.NORM_FROBENIUS) < 1e-14
+
+
 def test_constant():
     cg1 = FunctionSpace(UnitSquareMesh(5, 5), "CG", 1)
     f = assemble(interpolate(Constant(1.0), cg1))
@@ -660,3 +667,59 @@ def test_interpolate_composition(mode):
         res_adj = assemble(u5)
         res_adj2 = assemble(interpolate(TestFunction(V5), conj(TestFunction(V1)) * dx))
         assert np.allclose(res_adj.dat.data_ro, res_adj2.dat.data_ro)
+
+
+@pytest.mark.parallel([1, 3])
+def test_interpolate_form():
+    mesh = UnitSquareMesh(5, 5)
+    V3 = FunctionSpace(mesh, "CG", 3)
+    V2 = FunctionSpace(mesh, "CG", 2)
+    V1 = FunctionSpace(mesh, "CG", 1)
+
+    V3_trial = TrialFunction(V3)
+    V2_test = TestFunction(V2)
+    V1_test = TestFunction(V1)
+    V2_dual_trial = TrialFunction(V2.dual())
+
+    two_form = inner(V3_trial, V2_test) * dx  # V3 x V2 -> R, equiv V3 -> V2^*
+    interp = interpolate(V1_test, two_form)  # V3 x V1 -> R, equiv V3 -> V1^*
+    assert interp.arguments() == (V1_test, V3_trial)
+    res1 = assemble(interp)
+
+    I = interpolate(V1_test, V2_dual_trial)  # V2^* x V1 -> R, equiv V2^* -> V1^*
+    interp2 = action(I, two_form)  # V3 -> V1^*
+    assert interp2.arguments() == (V1_test, V3_trial)
+    res2 = assemble(interp2)
+    assert mat_equals(res1, res2)
+
+    res3 = assemble(inner(V3_trial, V1_test) * dx)  # V3 x V1 -> R
+    assert mat_equals(res1, res3)
+
+
+@pytest.mark.parallel([1, 3])
+def test_interpolate_form_mixed():
+    mesh = UnitSquareMesh(3, 3)
+    V1 = FunctionSpace(mesh, "CG", 1)
+    V2 = FunctionSpace(mesh, "CG", 2)
+    V3 = FunctionSpace(mesh, "CG", 3)
+    V4 = FunctionSpace(mesh, "CG", 4)
+    V = V3 * V4
+    W = V1 * V2
+
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    q = TestFunction(W)
+
+    form = inner(u, v) * dx  # V x V -> R, equiv V -> V^*
+    interp = interpolate(q, form)  # V -> W^*, equiv V x W -> R
+    assert interp.arguments() == (q, u)
+    res1 = assemble(interp)
+
+    I = interpolate(q, TrialFunction(V.dual()))  # V^* x W -> R, equiv V^* -> W^*
+    interp2 = action(I, form)  # V -> W^*
+    assert interp2.arguments() == (q, u)
+    res2 = assemble(interp2)
+    assert mat_equals(res1, res2)
+
+    res3 = assemble(inner(u, q) * dx)  # V x W -> R
+    assert mat_equals(res1, res3)