Move __future__ interpolation into interpolation.py

demos/benney_luke/benney_luke.py.rst

@@ -177,14 +177,8 @@ and define the exact solutions, which need to be updated at every time-step::
   expr_eta = 1/3.0*c*pow(cosh(0.5*sqrt(c*epsilon/mu)*(x[0]-x0-t_-epsilon*c*t_/6.0)),-2)
   expr_phi = 2/3.0*sqrt(c*mu/epsilon)*(tanh(0.5*sqrt(c*epsilon/mu)*(x[0]-x0-t_-epsilon*c*t_/6.0))+1)
 
-Since we will interpolate these values again and again, we use an
-:class:`~.Interpolator` whose :meth:`~.Interpolator.interpolate`
-method we can call to perform the interpolation. ::
-
-  eta_interpolator = Interpolator(expr_eta, ex_eta)
-  phi_interpolator = Interpolator(expr_phi, ex_phi)
-  phi_interpolator.interpolate()
-  eta_interpolator.interpolate()
+  ex_eta.interpolate(expr_eta)
+  ex_phi.interpolate(expr_phi)
 
 For visualisation, we save the computed and exact solutions to
 an output file.  Note that the visualised data will be interpolated
@@ -201,8 +195,8 @@ We are now ready to enter the main time iteration loop::
 
         t_.assign(t)
 
-        eta_interpolator.interpolate()
-        phi_interpolator.interpolate()
+        ex_eta.interpolate(expr_eta)
+        ex_phi.interpolate(expr_phi)
 
         phi_solver_h.solve()
         q_solver_h.solve()

demos/full_waveform_inversion/full_waveform_inversion.py.rst

@@ -217,8 +217,6 @@ The receiver mesh is required in order to interpolate the wave equation solution
 
 We are now able to proceed with the synthetic data computations and record them on the receivers::
 
-    from firedrake.__future__ import interpolate
-
     true_data_receivers = []
     total_steps = int(final_time / dt) + 1
     f = Cofunction(V.dual())  # Wave equation forcing term.

demos/higher_order_mass_lumping/higher_order_mass_lumping.py.rst

@@ -135,7 +135,7 @@ Using Firedrake, we specify the mass matrix using the special quadrature rule wi
     m = (u - 2.0 * u_n + u_nm1) / Constant(dt * dt) * v * dxlump
 
 .. note::
-    Mass lumping is a common technique in finite elements to produce a diagonal mass matrix that can be trivially inverted resulting in a in very efficient explicit time integration scheme. It's usually done with nodal basis functions and an inexact quadrature rule for the mass matrix. A diagonal matrix is obtained when the integration points coincide with the nodes of the basis function. However, when using elements of :math:`p \ge 2`, this technique does not result in a stable and accurate finite element scheme and new elements must be found such as those detailed in :cite:Chin:1999 .
+    Mass lumping is a common technique in finite elements to produce a diagonal mass matrix that can be trivially inverted resulting in a very efficient explicit time integration scheme. It's usually done with nodal basis functions and an inexact quadrature rule for the mass matrix. A diagonal matrix is obtained when the integration points coincide with the nodes of the basis function. However, when using elements of :math:`p \ge 2`, this technique does not result in a stable and accurate finite element scheme and new elements must be found such as those detailed in :cite:Chin:1999 .
 
 The stiffness matrix :math:`a(u,v)` is formed using a standard quadrature rule and is treated explicitly::
 

docs/notebooks/01-spd-helmholtz.ipynb

@@ -129,8 +129,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from firedrake import *\n",
-    "from firedrake.__future__ import interpolate"
+    "from firedrake import *"
    ]
   },
   {

docs/notebooks/03-elasticity.ipynb

@@ -187,7 +187,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "displaced_coordinates = interpolate(SpatialCoordinate(mesh) + uh, V)\n",
+    "displaced_coordinates = assemble(interpolate(SpatialCoordinate(mesh) + uh, V))\n",
     "displaced_mesh = Mesh(displaced_coordinates)"
    ]
   },

docs/notebooks/06-pde-constrained-optimisation.ipynb

@@ -226,7 +226,7 @@
     "u_inflow = as_vector([y*(10-y)/25.0, 0])\n",
     "\n",
     "noslip = DirichletBC(W.sub(0), (0, 0), (3, 5))\n",
-    "inflow = DirichletBC(W.sub(0), interpolate(u_inflow, V), 1)\n",
+    "inflow = DirichletBC(W.sub(0), assemble(interpolate(u_inflow, V)), 1)\n",
     "static_bcs = [inflow, noslip]"
    ]
   },

docs/notebooks/07-geometric-multigrid.ipynb

@@ -211,7 +211,7 @@
     "                       0)\n",
     "\n",
     "    value = as_vector([gbar*(1 - (12*t)**2), 0])\n",
-    "    bcs = [DirichletBC(W.sub(0), interpolate(value, V), (1, 2)),\n",
+    "    bcs = [DirichletBC(W.sub(0), assemble(interpolate(value, V)), (1, 2)),\n",
     "           DirichletBC(W.sub(0), zero(2), (3, 4))]\n",
     "    \n",
     "    a = (nu*inner(grad(u), grad(v)) - p*div(v) + q*div(u))*dx\n",
@@ -524,7 +524,7 @@
     "                       0)\n",
     "\n",
     "    value = as_vector([gbar*(1 - (12*t)**2), 0])\n",
-    "    bcs = [DirichletBC(W.sub(0), interpolate(value, V), (1, 2)),\n",
+    "    bcs = [DirichletBC(W.sub(0), assemble(interpolate(value, V)), (1, 2)),\n",
     "           DirichletBC(W.sub(0), zero(2), (3, 4))]\n",
     "    \n",
     "    a = (nu*inner(grad(u), grad(v)) - p*div(v) + q*div(u))*dx\n",

docs/notebooks/11-extract-adjoint-solutions.ipynb

@@ -133,7 +133,7 @@
     "\n",
     "# Assign initial condition\n",
     "x, y = SpatialCoordinate(mesh)\n",
-    "u_ = interpolate(as_vector([sin(pi*x), 0]), V)\n",
+    "u_ = assemble(interpolate(as_vector([sin(pi*x), 0]), V))\n",
     "u.assign(u_)\n",
     "\n",
     "# Set diffusivity constant\n",

docs/source/interpolation.rst

@@ -38,18 +38,6 @@ where :math:`\bar{\phi}^*_i` is the :math:`i`-th dual basis function to
 The interpolate operator
 ------------------------
 
-.. note::
-   The semantics for interpolation in Firedrake are in the course of changing.
-   The documentation provided here is for the new behaviour, in which the
-   `interpolate` operator is symbolic. In order to access the behaviour
-   documented here (which is recommended), users need to use the following
-   import line:
-
-   .. code-block:: python3
-
-      from firedrake.__future__ import interpolate
-
-
 The basic syntax for interpolation is:
 
 .. literalinclude:: ../../tests/firedrake/regression/test_interpolation_manual.py

firedrake/__future__.py

@@ -1,46 +1,16 @@
-from ufl.domain import as_domain, extract_unique_domain
-from ufl.algorithms import extract_arguments
-from firedrake.mesh import MeshTopology, VertexOnlyMeshTopology
-from firedrake.interpolation import (interpolate as interpolate_old,
-                                     Interpolator as InterpolatorOld,
-                                     SameMeshInterpolator as SameMeshInterpolatorOld,
-                                     CrossMeshInterpolator as CrossMeshInterpolatorOld)
-from firedrake.cofunction import Cofunction
-from functools import wraps
+from firedrake.interpolation import interpolate as interpolate_default, Interpolator as Interpolator_default
+import warnings
 
 
-__all__ = ("interpolate", "Interpolator")
-
-
-class Interpolator(InterpolatorOld):
-    def __new__(cls, expr, V, **kwargs):
-        target_mesh = as_domain(V)
-        source_mesh = extract_unique_domain(expr) or target_mesh
-        if target_mesh is source_mesh or all(isinstance(m.topology, MeshTopology) for m in [target_mesh, source_mesh]) and target_mesh.submesh_ancesters[-1] is source_mesh.submesh_ancesters[-1]:
-            return object.__new__(SameMeshInterpolator)
-        else:
-            if isinstance(target_mesh.topology, VertexOnlyMeshTopology):
-                return object.__new__(SameMeshInterpolator)
-            else:
-                return object.__new__(CrossMeshInterpolator)
-
-    interpolate = InterpolatorOld._interpolate_future
-
-
-class SameMeshInterpolator(Interpolator, SameMeshInterpolatorOld):
-    pass
-
-
-class CrossMeshInterpolator(Interpolator, CrossMeshInterpolatorOld):
-    pass
+def interpolate(expr, V, *args, **kwargs):
+    warnings.warn("""The symbolic `interpolate` has been moved into `firedrake`
+                  and is now the default method for interpolating in Firedrake. The need to
+                  `from firedrake.__future__ import interpolate` is now deprecated.""", FutureWarning)
+    return interpolate_default(expr, V, *args, **kwargs)
 
 
-@wraps(interpolate_old)
-def interpolate(expr, V, *args, **kwargs):
-    default_missing_val = kwargs.pop("default_missing_val", None)
-    if isinstance(V, Cofunction):
-        adjoint = bool(extract_arguments(expr))
-        return Interpolator(
-            expr, V.function_space().dual(), *args, **kwargs
-        ).interpolate(V, adjoint=adjoint, default_missing_val=default_missing_val)
-    return Interpolator(expr, V, *args, **kwargs).interpolate(default_missing_val=default_missing_val)
+class Interpolator(Interpolator_default):
+    def __new__(cls, *args, **kwargs):
+        warnings.warn("""The symbolic `Interpolator` has been moved into `firedrake`.
+                      The need to `from firedrake.__future__ import Interpolator` is now deprecated.""", FutureWarning)
+        return Interpolator_default(*args, **kwargs)

firedrake/interpolation.py

@@ -24,9 +24,10 @@
 import firedrake
 from firedrake import tsfc_interface, utils, functionspaceimpl
 from firedrake.ufl_expr import Argument, action, adjoint as expr_adjoint
-from firedrake.mesh import MissingPointsBehaviour, VertexOnlyMeshMissingPointsError
+from firedrake.mesh import MissingPointsBehaviour, VertexOnlyMeshMissingPointsError, VertexOnlyMeshTopology
 from firedrake.petsc import PETSc
 from firedrake.halo import _get_mtype as get_dat_mpi_type
+from firedrake.cofunction import Cofunction
 from mpi4py import MPI
 
 from pyadjoint import stop_annotating, no_annotations
@@ -142,20 +143,13 @@ def _ufl_expr_reconstruct_(self, expr, v=None, **interp_data):
 
 
 @PETSc.Log.EventDecorator()
-def interpolate(
-    expr,
-    V,
-    subset=None,
-    access=op2.WRITE,
-    allow_missing_dofs=False,
-    default_missing_val=None,
-    ad_block_tag=None
-):
-    """Interpolate an expression onto a new function in V.
+def interpolate(expr, V, *args, **kwargs):
+    """Returns a UFL expression for the interpolation operation of ``expr`` into ``V``.
 
     :arg expr: a UFL expression.
     :arg V: the :class:`.FunctionSpace` to interpolate into (or else
-        an existing :class:`.Function`).
+        an existing :class:`.Function` or :class:`.Cofunction`).
+        Adjoint interpolation requires ``V`` to be a :class:`.Cofunction`.
     :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`.
@@ -182,8 +176,7 @@ def interpolate(
         to zero. Ignored if interpolating within the same mesh or onto a
         :func:`.VertexOnlyMesh`.
     :kwarg ad_block_tag: An optional string for tagging the resulting assemble block on the Pyadjoint tape.
-    :returns: a new :class:`.Function` in the space ``V`` (or ``V`` if
-        it was a Function).
+    :returns: A synbolic :class:`.Interpolate` object
 
     .. note::
 
@@ -204,9 +197,13 @@ def interpolate(
        performance by using an :class:`Interpolator` instead.
 
     """
-    return Interpolator(
-        expr, V, subset=subset, access=access, allow_missing_dofs=allow_missing_dofs
-    ).interpolate(default_missing_val=default_missing_val, ad_block_tag=ad_block_tag)
+    default_missing_val = kwargs.pop("default_missing_val", None)
+    if isinstance(V, Cofunction):
+        adjoint = bool(extract_arguments(expr))
+        return Interpolator(
+            expr, V.function_space().dual(), *args, **kwargs
+        ).interpolate(V, adjoint=adjoint, default_missing_val=default_missing_val)
+    return Interpolator(expr, V, *args, **kwargs).interpolate(default_missing_val=default_missing_val)
 
 
 class Interpolator(abc.ABC):
@@ -264,7 +261,7 @@ def __new__(cls, expr, V, **kwargs):
         if target_mesh is source_mesh or submesh_interp_implemented:
             return object.__new__(SameMeshInterpolator)
         else:
-            if isinstance(target_mesh.topology, firedrake.mesh.VertexOnlyMeshTopology):
+            if isinstance(target_mesh.topology, VertexOnlyMeshTopology):
                 return object.__new__(SameMeshInterpolator)
             else:
                 return object.__new__(CrossMeshInterpolator)
@@ -296,7 +293,7 @@ def __init__(
             expr = replace(expr, {v: v.reconstruct(number=1)})
         self.expr_renumbered = expr
 
-    def _interpolate_future(self, *function, transpose=None, adjoint=False, default_missing_val=None):
+    def interpolate(self, *function, transpose=None, adjoint=False, default_missing_val=None):
         """Define the :class:`Interpolate` object corresponding to the interpolation operation of interest.
 
         Parameters
@@ -322,11 +319,6 @@ def _interpolate_future(self, *function, transpose=None, adjoint=False, default_
         -------
         firedrake.interpolation.Interpolate or ufl.action.Action or ufl.adjoint.Adjoint
             The symbolic object representing the interpolation operation.
-
-        Notes
-        -----
-        This method is the default future behaviour of interpolation. In a future release, the
-        ``Interpolator.interpolate`` method will be replaced by this method.
         """
 
         V = self.V
@@ -351,79 +343,6 @@ def _interpolate_future(self, *function, transpose=None, adjoint=False, default_
         # Return the `ufl.Interpolate` object
         return interp
 
-    @PETSc.Log.EventDecorator()
-    def interpolate(self, *function, output=None, transpose=None, adjoint=False, default_missing_val=None,
-                    ad_block_tag=None):
-        """Compute the interpolation by assembling the appropriate :class:`Interpolate` object.
-
-        Parameters
-        ----------
-        *function: firedrake.function.Function or firedrake.cofunction.Cofunction
-                   If the expression being interpolated contains an argument,
-                   then the function value to interpolate.
-        output: firedrake.function.Function or firedrake.cofunction.Cofunction
-                A function to contain the output.
-        transpose : bool
-                   Deprecated, use adjoint instead.
-        adjoint: bool
-                   Set to true to apply the adjoint of the interpolation
-                   operator.
-        default_missing_val: bool
-                             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 the values are
-                             either (a) unchanged if some ``output`` is specified to the
-                             :meth:`interpolate` method or (b) set to zero. This does not affect
-                             adjoint interpolation. Ignored if interpolating within the same
-                             mesh or onto a :func:`.VertexOnlyMesh`.
-        ad_block_tag: str
-                      An optional string for tagging the resulting assemble block on the Pyadjoint tape.
-
-        Returns
-        -------
-        firedrake.function.Function or firedrake.cofunction.Cofunction
-            The resulting interpolated function.
-        """
-        from firedrake.assemble import assemble
-
-        warnings.warn("""The use of `interpolate` to perform the numerical interpolation is deprecated.
-This feature will be removed very shortly.
-
-Instead, import `interpolate` from the `firedrake.__future__` module to update
-the interpolation's behaviour to return the symbolic `ufl.Interpolate` object associated
-with this interpolation.
-
-You can then assemble the resulting object to get the interpolated quantity
-of interest. For example,
-
-```
-from firedrake.__future__ import interpolate
-...
-
-assemble(interpolate(expr, V))
-```
-
-Alternatively, you can also perform other symbolic operations on the interpolation operator, such as taking
-the derivative, and then assemble the resulting form.
-""", FutureWarning)
-        if transpose is not None:
-            warnings.warn("'transpose' argument is deprecated, use 'adjoint' instead", FutureWarning)
-            adjoint = transpose or adjoint
-
-        # Get the Interpolate object
-        interp = self._interpolate_future(*function, adjoint=adjoint,
-                                          default_missing_val=default_missing_val)
-
-        if isinstance(self.V, firedrake.Function) and not output:
-            # V can be the Function to interpolate into (e.g. see `Function.interpolate``).
-            output = self.V
-
-        # Assemble the `ufl.Interpolate` object, which will then call `Interpolator._interpolate`
-        # to perform the interpolation. Having this structure ensures consistency between
-        # `Interpolator` and `Interp`. This mechanism handles annotation since performing interpolation will drop an
-        # `AssembleBlock` on the tape.
-        return assemble(interp, tensor=output, ad_block_tag=ad_block_tag)
-
     @abc.abstractmethod
     def _interpolate(self, *args, **kwargs):
         """
@@ -715,10 +634,10 @@ def _interpolate(
                     # so the sub_interpolators are already prepared to interpolate
                     # without needing to be given a Function
                     assert not self.nargs
-                    interp = sub_interpolator._interpolate_future(adjoint=adjoint, **kwargs)
+                    interp = sub_interpolator.interpolate(adjoint=adjoint, **kwargs)
                     assemble(interp, tensor=output_sub_func)
                 else:
-                    interp = sub_interpolator._interpolate_future(adjoint=adjoint, **kwargs)
+                    interp = sub_interpolator.interpolate(adjoint=adjoint, **kwargs)
                     assemble(action(interp, f_src_sub_func), tensor=output_sub_func)
             return output
 

firedrake/preconditioners/hypre_ads.py

@@ -5,7 +5,7 @@
 from firedrake.ufl_expr import TestFunction
 from firedrake.dmhooks import get_function_space
 from firedrake.preconditioners.hypre_ams import chop
-from firedrake.__future__ import interpolate
+from firedrake.interpolation import interpolate
 from ufl import grad, curl, SpatialCoordinate
 from pyop2.utils import as_tuple