Interpolate: Fix constant and coefficient ordering

Author: pbrubeck

firedrake/assemble.py

@@ -553,8 +553,7 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
             arg_expression = ufl.algorithms.extract_arguments(expression)
 
             # Handle interpolation onto subfunctions
-            if (arg_expression and len(arg_expression[0].function_space()) > 1
-                    and not isinstance(expression, firedrake.Argument)):
+            if arg_expression and len(arg_expression[0].function_space()) > 1:
                 assert rank == 1
                 V = arg_expression[0].function_space()
                 if tensor is not None:

firedrake/interpolation.py

@@ -1196,6 +1196,10 @@ def _interpolator(V, tensor, expr, subset, arguments, access, bcs=None):
     if needs_cell_sizes:
         cs = target_mesh.cell_sizes
         parloop_args.append(cs.dat(op2.READ, cs.cell_node_map()))
+
+    for const in extract_firedrake_constants(expr):
+        parloop_args.append(const.dat(op2.READ))
+
     for coefficient in coefficients:
         if isinstance(target_mesh.topology, firedrake.mesh.VertexOnlyMeshTopology):
             coeff_mesh = extract_unique_domain(coefficient)
@@ -1224,9 +1228,6 @@ def _interpolator(V, tensor, expr, subset, arguments, access, bcs=None):
             m_ = coefficient.function_space().entity_node_map(target_mesh.topology, "cell", None, None)
         parloop_args.append(coefficient.dat(op2.READ, m_))
 
-    for const in extract_firedrake_constants(expr):
-        parloop_args.append(const.dat(op2.READ))
-
     parloop = op2.ParLoop(*parloop_args)
     parloop_compute_callable = parloop.compute
     if isinstance(tensor, op2.Mat):

tests/firedrake/regression/test_interp_dual.py

@@ -280,7 +280,7 @@ def test_solve_interp_u(mesh):
 
 
 @pytest.mark.parametrize("family0,degree0,family1,degree1",
-                         [("DG", 1, "CG", 1),
+                         [("DG", 1, "CG", 2),
                           ("DG", 0, "RT", 1)])
 def test_interp_subfunction(mesh, family0, degree0, family1, degree1):
     V = FunctionSpace(mesh, "DG", 0)
@@ -295,7 +295,7 @@ def test_interp_subfunction(mesh, family0, degree0, family1, degree1):
 
     expr = sum(w[i] for i in np.ndindex(w.ufl_shape))
 
-    Fw = inner(1, expr)*dx
+    Fw = inner(1, expr)*dx(degree=0)
     expected = assemble(Fw)
 
     IFv = Interpolate(expr, Fv)
@@ -306,5 +306,7 @@ def test_interp_subfunction(mesh, family0, degree0, family1, degree1):
 
         result = assemble(IFv, tensor=tensor)
         assert result.function_space() == W.dual()
-        with result.dat.vec_ro as x, expected.dat.vec_ro as y:
-            assert np.allclose(x[:], y[:])
+        if tensor:
+            assert result is tensor
+        for x, y, in zip(result.subfunctions, expected.subfunctions):
+            assert np.allclose(x.dat.data_ro, y.dat.data_ro)