add interpolation tests

firedrake/interpolation.py

@@ -905,6 +905,8 @@ def make_interpolator(expr, V, subset, access, bcs=None):
     elif len(arguments) == 1:
         if isinstance(V, firedrake.Function):
             raise ValueError("Cannot interpolate an expression with an argument into a Function")
+        if len(V) > 1:
+            raise NotImplementedError("Interpolation of mixed expressions with arguments is not supported")
         argfs = arguments[0].function_space()
         source_mesh = argfs.mesh()
         argfs_map = argfs.cell_node_map()
@@ -995,23 +997,25 @@ def callable():
 
         if len(V) == 1:
             loops.extend(_interpolator(V, tensor, expr, subset, arguments, access, bcs=bcs))
-        elif (hasattr(expr, "subfunctions") and len(expr.subfunctions) == len(V)
-              and all(sub_expr.ufl_shape == Vsub.value_shape for Vsub, sub_expr in zip(V, expr.subfunctions))):
-            for Vsub, sub_tensor, sub_expr in zip(V, tensor, expr.subfunctions):
-                loops.extend(_interpolator(Vsub, sub_tensor, sub_expr, subset, arguments, access, bcs=bcs))
-        elif len(arguments) == 0:
-            # Unflatten the expression into each of the mixed components
-            offset = 0
-            for Vsub, sub_tensor in zip(V, tensor):
-                if len(Vsub.value_shape) == 0:
-                    sub_expr = expr[offset]
-                else:
-                    components = [expr[offset + j] for j in range(Vsub.value_size)]
-                    sub_expr = ufl.as_tensor(numpy.reshape(components, Vsub.value_shape))
-                loops.extend(_interpolator(Vsub, sub_tensor, sub_expr, subset, arguments, access, bcs=bcs))
-                offset += Vsub.value_size
         else:
-            raise NotImplementedError("Cannot interpolate a mixed expression with %d arguments" % len(arguments))
+            if (hasattr(expr, "subfunctions") and len(expr.subfunctions) == len(V)
+                    and all(sub_expr.ufl_shape == Vsub.value_shape for Vsub, sub_expr in zip(V, expr.subfunctions))):
+                # Use subfunctions if they match the target shapes
+                expressions = expr.subfunctions
+            else:
+                # Unflatten the expression into the shapes of the mixed components
+                offset = 0
+                expressions = []
+                for Vsub in V:
+                    if len(Vsub.value_shape) == 0:
+                        expressions.append(expr[offset])
+                    else:
+                        components = [expr[offset + j] for j in range(Vsub.value_size)]
+                        expressions.append(ufl.as_tensor(numpy.reshape(components, Vsub.value_shape)))
+                    offset += Vsub.value_size
+            # Interpolate each sub expression into each function space
+            for Vsub, sub_tensor, sub_expr in zip(V, tensor, expressions):
+                loops.extend(_interpolator(Vsub, sub_tensor, sub_expr, subset, arguments, access, bcs=bcs))
 
         if bcs and len(arguments) == 0:
             loops.extend(partial(bc.apply, f) for bc in bcs)

tests/regression/test_interpolate.py

@@ -28,6 +28,47 @@ def test_function():
     assert np.allclose(g.dat.data, h.dat.data)
 
 
+def test_mixed_expression():
+    m = UnitTriangleMesh()
+    x = SpatialCoordinate(m)
+    V1 = FunctionSpace(m, 'P', 1)
+    V2 = FunctionSpace(m, 'P', 2)
+
+    V = V1 * V2
+    expressions = [x[0], x[0]*x[1]]
+    expr = as_vector(expressions)
+    fg = assemble(interpolate(expr, V))
+    f, g = fg.subfunctions
+
+    f1 = Function(V1).interpolate(expressions[0])
+    g1 = Function(V2).interpolate(expressions[1])
+    assert np.allclose(f.dat.data, f1.dat.data)
+    assert np.allclose(g.dat.data, g1.dat.data)
+
+
+def test_mixed_function():
+    m = UnitTriangleMesh()
+    x = SpatialCoordinate(m)
+    V1 = FunctionSpace(m, 'RT', 1)
+    V2 = FunctionSpace(m, 'DG', 0)
+    V = V1 * V2
+
+    expressions = [x[0], x[1], Constant(0.444)]
+    expr = as_vector(expressions)
+    v = assemble(interpolate(expr, V))
+
+    W1 = FunctionSpace(m, 'RT', 2)
+    W2 = FunctionSpace(m, 'DG', 1)
+    W = W1 * W2
+    w = assemble(interpolate(v, W))
+
+    f, g = w.subfunctions
+    f1 = Function(W1).interpolate(x)
+    g1 = Function(W2).interpolate(expressions[-1])
+    assert np.allclose(f.dat.data, f1.dat.data)
+    assert np.allclose(g.dat.data, g1.dat.data)
+
+
 def test_inner():
     m = UnitTriangleMesh()
     V1 = FunctionSpace(m, 'P', 1)

tests/vertexonly/test_interpolation_from_parent.py

@@ -292,13 +292,19 @@ def test_tensor_function_interpolation(parentmesh, vertexcoords, tfs):
     assert np.allclose(w_v.dat.data_ro.reshape(result.shape), 2*result)
 
 
-@pytest.mark.xfail(raises=NotImplementedError, reason="Interpolation of UFL expressions into mixed functions not supported")
 def test_mixed_function_interpolation(parentmesh, vertexcoords, tfs):
     if parentmesh.name == "immersedsphere":
         vertexcoords = immersed_sphere_vertexcoords(parentmesh, vertexcoords)
     tfs_fam, tfs_deg, tfs_typ = tfs
+
     vm = VertexOnlyMesh(parentmesh, vertexcoords, missing_points_behaviour=None)
-    vertexcoords = vm.coordinates.dat.data_ro.reshape(-1, parentmesh.geometric_dimension())
+    vertexcoords = vm.coordinates.dat.data_ro
+    if (
+        parentmesh.coordinates.function_space().ufl_element().family()
+        == "Discontinuous Lagrange"
+    ):
+        pytest.skip(f"Interpolating into f{tfs_fam} on a periodic mesh is not well-defined.")
+
     V1 = tfs_typ(parentmesh, tfs_fam, tfs_deg)
     V2 = FunctionSpace(parentmesh, "CG", 1)
     V = V1 * V2
@@ -311,30 +317,23 @@ def test_mixed_function_interpolation(parentmesh, vertexcoords, tfs):
     # Get Function in V1
     # use outer product to check Regge works
     expr1 = outer(x, x)
-    assert W1.shape == expr1.ufl_shape
-    assemble(interpolate(expr1, v1))
+    assert W1.value_shape == expr1.ufl_shape
+    v1.interpolate(expr1)
     result1 = np.asarray([np.outer(vertexcoords[i], vertexcoords[i]) for i in range(len(vertexcoords))])
     if len(result1) == 0:
         result1 = result1.reshape(vertexcoords.shape + (parentmesh.geometric_dimension(),))
     # Get Function in V2
     expr2 = reduce(add, SpatialCoordinate(parentmesh))
-    assemble(interpolate(expr2, v2))
+    v2.interpolate(expr2)
     result2 = np.sum(vertexcoords, axis=1)
+
     # Interpolate Function in V into W
     w_v = assemble(interpolate(v, W))
+
     # Split result and check
     w_v1, w_v2 = w_v.subfunctions
-    assert np.allclose(w_v1.dat.data_ro, result1)
-    assert np.allclose(w_v2.dat.data_ro, result2)
-    # try and make reusable Interpolator from V to W
-    A_w = Interpolator(TestFunction(V), W)
-    w_v = Function(W)
-    assemble(A_w.interpolate(v), tensor=w_v)
-    # Split result and check
-    w_v1, w_v2 = w_v.subfunctions
-    assert np.allclose(w_v1.dat.data_ro, result1)
-    assert np.allclose(w_v2.dat.data_ro, result2)
-    # Enough tests - don't both using it again for a different Function in V
+    assert np.allclose(w_v1.dat.data_ro.reshape(result1.shape), result1)
+    assert np.allclose(w_v2.dat.data_ro.reshape(result2.shape), result2)
 
 
 def test_scalar_real_interpolation(parentmesh, vertexcoords):
@@ -431,7 +430,6 @@ def test_tensor_function_interpolation_parallel(parentmesh, vertexcoords, tfs):
     test_tensor_function_interpolation(parentmesh, vertexcoords, tfs)
 
 
-@pytest.mark.xfail(reason="Interpolation of UFL expressions into mixed functions not supported")
 @pytest.mark.parallel
 def test_mixed_function_interpolation_parallel(parentmesh, vertexcoords, tfs):
     test_mixed_function_interpolation(parentmesh, vertexcoords, tfs)