Implement cross mesh interpolation

This implements cross-mesh interpolation for interpolating into function
spaces which have point evaluation nodes. Full documentation is added to
the manual.

Two new classes are added to interpolation.py: CrossMeshInterpolator and
SameMeshInterpolator, whilst Interpolator is made an abstract base
class. To maintain the API of interpolate and Interpolator, the __new__
method of Interpolator is  overridden to return an instance of the
appropriate subclass.

Two new keyword arguments are added to interpolate and Interpolator to
allow for target meshes which extend outside the source mesh domain::
see their docstrings for details.

Docstrings are also added for some undocumented keyword arguments of
Interpolator and interpolate.

A full suite of tests is found in
tests/regression/test_interpolate_cross_mesh.py

Note that as part of this, I have changed the error by VertexOnlyMesh
when points are outside the domain from a ValueError to a
VertexOnlyMeshMissingPointsError

Details of how this works from the PR description:

We use VertexOnlyMesh as an intermediary for the global locations of the
point evaluation nodes of the target function space:

1. We get the point evaluation node locations using the method described
in the interpolation from external data section of the manual. This will
have the parallel domain decomposition of the target mesh.
2. Next we create a VertexOnlyMesh (A) at those locations within the
source mesh such that we inherit the source mesh's parallel domain
decomposition.
3. We interpolate our expression in our source function space onto a
P0DG function space on VertexOnlyMesh (A), which has the effect of point
evaluating at the target function space node locations.
4. This VertexOnlyMesh (A) has an input_ordering VertexOnlyMesh (B)
whose vertices have the ordering and parallel domain decomposition of
the target function space global node locations. We interpolate from
P0DG on (A) onto P0DG on (B). Under the hood, this is an SF reduce
operation which moves the point evaluations from (A) to (B).
5. We now have a Function on the input_ordering VertexOnlyMesh (B) which
has point evaluations from our source mesh function space at the target
mesh function space node locations. These are in the correct order and
have the correct domain decomposition. We can therefore safely set the
dat of a function in our target function space to the values of this
function.

For this to work for the general case, we would need to create a
VertexOnlyMesh at the global quadrature points of the target function
space, which is rather more complicated than the work I've done here.

Some important notes:

 - This does not require one mesh to be a structured refinement of the
 other. This, for example, should allow you to solve a PDE on two
 different unstructured meshes, one of which is finer than the other,
 and directly check the difference in solutions by interpolating from
 one mesh to the other within Firedrake.
 - Crucially, this is entirely parallel compatible!
 - Since we can interpolate onto immersed manifolds we can perform line,
 surface and volume integrals by interpolating onto a mesh which has the
 domain of integration as an immersed manifold. This is demonstrated in
 the test_interpolate_cross_mesh.py test suite.

Other notes:

- The VertexOnlyMesh required is stored in the Interpolator, as are the
 underlying Interpolators
- For interpolation between mixed spaces, I create sub_interpolators for
each space and evaluate them as necessary when calling interpolate
- Interpolation into a mixed space therefore requires the function space
 being interpolated from to be another mixed space. I throw a
 NotImplementedError if not.

Regarding the manual:
Note that not all the comments in the manual file are included in the
literalinclude text of the manual, instead they are approximately
rewritten as prose in the manual.

docs/source/interpolation.rst

@@ -108,6 +108,120 @@ in the dual space to `V` to assembled 1-forms in the dual space to
 `W`.
 
 
+Interpolation across meshes
+---------------------------
+
+The interpolation API supports interpolation between meshes where the target
+function space has finite elements (as given in the list of
+:ref:`supported elements <supported_elements>`)
+
+* **Lagrange/CG** (also known a Continuous Galerkin or P elements),
+* **Q** (i.e. Lagrange/CG on lines, quadrilaterals and hexahedra),
+* **Discontinuous Lagrange/DG** (also known as Discontinuous Galerkin or DP elements) and
+* **DQ** (i.e. Discontinuous Lagrange/DG on lines, quadrilaterals and hexahedra).
+
+Vector, tensor and mixed function spaces can also be interpolated into from
+other meshes as long as they are constructed from these spaces.
+
+.. note::
+
+   The list of supported elements above is only for *target* function spaces.
+   Function spaces on the *source* mesh can be built from most of the supported
+   elements.
+
+There are few constraints on the meshes involved: the target mesh can have a
+different cell shape, topological dimension, or resolution to the source mesh.
+There are many use cases for this: For example, two solutions to the same
+problem calculated on meshes with different resolutions or cell shapes can be
+interpolated onto one another, or onto a third, finer mesh, and be directly
+compared.
+
+
+Interpolating onto sub-domain meshes
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The target mesh for a cross-mesh interpolation need not cover the full domain
+of the source mesh. Volume, surface and line integrals can therefore be
+calculated by interpolating onto the mesh or
+:ref:`immersed manifold <immersed_manifolds>` which defines the volume,
+surface or line of interest in the domain. The integral itself is calculated
+by calling :py:func:`~.assemble` on an approriate form over the target mesh
+function space:
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :lines: 11-18, 31-40
+
+For more on forms, see :ref:`this section of the manual <more_complicated_forms>`.
+
+
+Interpolating onto other meshes
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. note::
+
+   Interpolation *from* :ref:`immersed manifolds <immersed_manifolds>` and
+   :ref:`high-order meshes <changing_coordinate_fs>` is currently not supported.
+
+If the target mesh extends outside the source mesh domain, then cross-mesh
+interpolation will raise a :py:class:`~.DofNotDefinedError`.
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :lines: 46-58, 64-65
+
+This can be overriden with the optional ``allow_missing_dofs`` keyword
+argument:
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :lines: 72-73, 80
+
+By default the missing degrees of freedom (DoFs, the global basis function
+coefficients which could not be set) are zero:
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :lines: 85
+
+We can optionally specify a value to use for our missing DoFs. Here
+we set them to be ``nan`` ('not a number') for easy identification:
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :lines: 90-93
+
+When using :py:class:`~.Interpolator`\s, the ``allow_missing_dofs`` keyword
+argument is set at construction:
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :lines: 100
+
+The ``default_missing_val`` keyword argument is then set whenever we call
+:py:meth:`~.Interpolator.interpolate`:
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :lines: 103
+
+If we supply an output :py:class:`~.Function` and don't set
+``default_missing_val`` then any missing DoFs are left as they were prior to
+interpolation:
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :lines: 107-110, 113-115, 124-126
+
+
 Interpolation from external data
 --------------------------------
 

docs/source/mesh-coordinates.rst

@@ -36,6 +36,8 @@ This can also be used if `f` is a solution to a PDE.
    This will be fixed in a future Firedrake update.
 
 
+.. _changing_coordinate_fs:
+
 Changing the coordinate function space
 --------------------------------------
 

docs/source/point-evaluation.rst

@@ -217,9 +217,10 @@ duplication.
 Points outside the domain
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 
-Be default points outside the domain  by more than the :ref:`specified
-tolerance <tolerance>` will generate a ``ValueError``. This can be switched
-to a warning or switched off entirely:
+Be default points outside the domain by more than the :ref:`specified
+tolerance <tolerance>` will generate
+a :class:`~.VertexOnlyMeshMissingPointsError`. This can be switched to a
+warning or switched off entirely:
 
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3

docs/source/variational-problems.rst

@@ -77,6 +77,7 @@ of standard shapes.  1-dimensional intervals may be constructed with
 example to build unit square meshes).  See
 :mod:`~firedrake.utility_meshes` for full details.
 
+.. _immersed_manifolds:
 
 Immersed manifolds
 ~~~~~~~~~~~~~~~~~~
@@ -260,6 +261,9 @@ the horizontal and vertical spaces when building a function space on
 an extruded mesh.  We refer the reader to the :doc:`manual section on
 extrusion <extruded-meshes>` for details.
 
+
+.. _supported_elements:
+
 Supported finite elements
 -------------------------
 
@@ -632,6 +636,8 @@ we can write:
        t += dt
        c.assign(t)
 
+.. _more_complicated_forms:
+
 More complicated forms
 ----------------------
 

firedrake/function.py

@@ -367,14 +367,47 @@ def vector(self):
         return vector.Vector(self)
 
     @PETSc.Log.EventDecorator()
-    def interpolate(self, expression, subset=None, ad_block_tag=None):
+    def interpolate(
+        self,
+        expression,
+        subset=None,
+        allow_missing_dofs=False,
+        default_missing_val=None,
+        ad_block_tag=None,
+    ):
         r"""Interpolate an expression onto this :class:`Function`.
 
         :param expression: a UFL expression to interpolate
-        :param ad_block_tag: string for tagging the resulting block on the Pyadjoint tape
+        :kwarg subset: An optional :class:`pyop2.types.set.Subset` to apply the
+            interpolation over. Cannot, at present, be used when interpolating
+            across meshes unless the target mesh is a :func:`.VertexOnlyMesh`.
+        :kwarg allow_missing_dofs: For interpolation across meshes: allow
+            degrees of freedom (aka DoFs/nodes) in the target mesh that cannot be
+            defined on the source mesh. For example, where nodes are point
+            evaluations, points in the target mesh that are not in the source mesh.
+            When ``False`` this raises a ``ValueError`` should this occur. When
+            ``True`` the corresponding values are set to zero or to the value
+            ``default_missing_val`` if given. Ignored if interpolating within the
+            same mesh or onto a :func:`.VertexOnlyMesh` (the behaviour of a
+            :func:`.VertexOnlyMesh` in this scenario is, at present, set when
+            it is created).
+        :kwarg default_missing_val: For interpolation across meshes: the optional
+            value to assign to DoFs in the target mesh that are outside the source
+            mesh. If this is not set then zero is used. Ignored if interpolating
+            within the same mesh or onto a :func:`.VertexOnlyMesh`.
+        :kwarg ad_block_tag: An optional string for tagging the resulting block on
+            the Pyadjoint tape.
         :returns: this :class:`Function` object"""
         from firedrake import interpolation
-        return interpolation.interpolate(expression, self, subset=subset, ad_block_tag=ad_block_tag)
+
+        return interpolation.interpolate(
+            expression,
+            self,
+            subset=subset,
+            allow_missing_dofs=allow_missing_dofs,
+            default_missing_val=default_missing_val,
+            ad_block_tag=ad_block_tag,
+        )
 
     def zero(self, subset=None):
         """Set all values to zero.