Merge pull request #3067 from firedrakeproject/ReubenHill/interpolation-mesh-to-mesh

Parallel Compatible Cross-Mesh Interpolation

Author: ReubenHill

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.