add docs

Author: leo-collins

docs/source/interpolation.rst

@@ -109,7 +109,7 @@ This also works when interpolating into a space defined on the facets of the mes
 
 
 Semantics of symbolic interpolation
---------------------------------------------
+-----------------------------------
 
 Let :math:`U` and :math:`V` be finite element spaces with DoFs :math:`\{\psi^{*}_{i}\}` and :math:`\{\phi^{*}_{i}\}`
 and basis functions :math:`\{\psi_{i}\}` and :math:`\{\phi_{i}\}`, respectively.
@@ -430,6 +430,63 @@ the next operation using ``f``.
 For interaction with external point data, see the
 :ref:`corresponding manual section <external-point-data>`.
 
+Interpolation between mixed function spaces
+-------------------------------------------
+
+Assembly of interpolation operators between mixed function spaces is also supported.
+Each component of the mixed space may be on different meshes.
+For example, consider the following mixed finite element spaces:
+
+.. math::
+
+   W &= V_1 \times V_2 \\
+   U &= V_3 \times V_4
+
+where each :math:`V_i` is a finite element space defined on possibly different meshes.
+We can assemble the interpolation matrix from :math:`U` to :math:`W` in Firedrake as follows:
+
+.. literalinclude:: ../../tests/firedrake/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_mixed_space_interpolation 1]
+   :end-before: [test_mixed_space_interpolation 2]
+
+We specified ``mat_type="nest"`` here to obtain a PETSc MatNest matrix, but Firedrake also
+supports assembly of ``mat_type="aij"`` and ``mat_type="matfree"`` interpolation matrices
+between mixed function spaces. In this example ``I`` is a block diagonal matrix, with
+each block given by
+
+.. math::
+
+   \begin{pmatrix}
+   V_3 \rightarrow V_1 & 0 \\
+   0 & V_4 \rightarrow V_2
+   \end{pmatrix}
+
+The off-diagonal blocks are zero since the dofs are applied component-wise. Firedrake's form
+compiler recognises this and avoids assembling the zero blocks.
+
+We can assemble more general interpolation matrices between mixed function spaces by interpolating
+vector expressions with arguments. For example, by doing
+
+.. literalinclude:: ../../tests/firedrake/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_mixed_space_interpolation 3]
+   :end-before: [test_mixed_space_interpolation 4]
+
+we can assemble the interpolation matrix with block structure
+
+.. math::
+
+   \begin{pmatrix}
+   V_3 \rightarrow V_1 & V_4 \rightarrow V_1 \\
+   V_3 \rightarrow V_2 & V_4 \rightarrow V_2
+   \end{pmatrix}
+
+Here we obtain non-zero off-diagonal blocks by including both components of the trial function
+in each component of the expression.
+
 Generating Functions with randomised values
 -------------------------------------------
 

tests/firedrake/regression/test_interpolation_manual.py

@@ -1,4 +1,5 @@
 from firedrake import *
+from firedrake.formmanipulation import split_form
 import pytest
 import numpy as np
 
@@ -235,3 +236,50 @@ def correct_indent():
 
     assert np.isclose(dest_eval05.evaluate(f_dest), 0.5)  # x_src^2 + y_src^2 = 0.5
     assert np.isclose(dest_eval15.evaluate(f_dest), 3.0)  # x_dest + y_dest = 3.0
+
+
+def test_mixed_space_interpolation():
+    mesh = UnitSquareMesh(2, 2)
+    V1 = FunctionSpace(mesh, "CG", 1)
+    V2 = FunctionSpace(mesh, "CG", 2)
+    V3 = FunctionSpace(mesh, "CG", 3)
+    V4 = FunctionSpace(mesh, "CG", 4)
+    W = V1 * V2
+    U = V3 * V4
+
+    # [test_mixed_space_interpolation 1]
+    interp = interpolate(TrialFunction(U), W)
+    I = assemble(interp, mat_type="nest")
+    # [test_mixed_space_interpolation 2]
+
+    # The block matrix structure is
+    # | V3 -> V1      0    |
+    # |     0     V4 -> V2 |
+    for i in range(2):
+        for j in range(2):
+            sub_mat = I.petscmat.getNestSubMatrix(i, j)
+            if i != j:
+                assert not sub_mat
+                continue
+            else:
+                res_block = assemble(interpolate(TrialFunction(U.sub(j)), W.sub(i)))
+                assert np.allclose(sub_mat[:, :], res_block.petscmat[:, :])
+                assert sub_mat.type == "seqaij"
+
+    # [test_mixed_space_interpolation 3]
+    u0, u1 = TrialFunctions(U)
+    expr = as_vector([u0 + u1, u0 + u1])
+    interp = interpolate(expr, W)
+    I2 = assemble(interp, mat_type="nest")
+    # [test_mixed_space_interpolation 4]
+
+    # The block matrix structure is
+    # | V3 -> V1   V4 -> V1 |
+    # | V3 -> V2   V4 -> V2 |
+    split_interp = dict(split_form(interp))
+    for i in range(2):
+        for j in range(2):
+            interp_ij = split_interp[(i, j)]
+            assert isinstance(interp_ij, Interpolate)
+            res_block = assemble(interpolate(TrialFunction(U.sub(j)), W.sub(i), allow_missing_dofs=True))
+            assert np.allclose(I2.petscmat.getNestSubMatrix(i, j)[:, :], res_block.petscmat[:, :])