diff --git a/demos/benney_luke/benney_luke.py.rst b/demos/benney_luke/benney_luke.py.rst
index 2931926249..0d035af52e 100644
--- a/demos/benney_luke/benney_luke.py.rst
+++ b/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()
diff --git a/demos/full_waveform_inversion/full_waveform_inversion.py.rst b/demos/full_waveform_inversion/full_waveform_inversion.py.rst
index eb4c135756..fcc2ffaadc 100644
--- a/demos/full_waveform_inversion/full_waveform_inversion.py.rst
+++ b/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.
diff --git a/demos/higher_order_mass_lumping/higher_order_mass_lumping.py.rst b/demos/higher_order_mass_lumping/higher_order_mass_lumping.py.rst
index 3f4e3bde28..193047aa27 100644
--- a/demos/higher_order_mass_lumping/higher_order_mass_lumping.py.rst
+++ b/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::
 
diff --git a/docs/notebooks/01-spd-helmholtz.ipynb b/docs/notebooks/01-spd-helmholtz.ipynb
index 021ce7fcfc..adb2849b60 100644
--- a/docs/notebooks/01-spd-helmholtz.ipynb
+++ b/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 *"
    ]
   },
   {
diff --git a/docs/notebooks/03-elasticity.ipynb b/docs/notebooks/03-elasticity.ipynb
index 15feea6d97..97d290d6f0 100644
--- a/docs/notebooks/03-elasticity.ipynb
+++ b/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)"
    ]
   },
diff --git a/docs/notebooks/06-pde-constrained-optimisation.ipynb b/docs/notebooks/06-pde-constrained-optimisation.ipynb
index 1b3cbda968..0a4e965448 100644
--- a/docs/notebooks/06-pde-constrained-optimisation.ipynb
+++ b/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]"
    ]
   },
diff --git a/docs/notebooks/07-geometric-multigrid.ipynb b/docs/notebooks/07-geometric-multigrid.ipynb
index 3eefe36deb..ace7b7caa2 100644
--- a/docs/notebooks/07-geometric-multigrid.ipynb
+++ b/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",
diff --git a/docs/notebooks/11-extract-adjoint-solutions.ipynb b/docs/notebooks/11-extract-adjoint-solutions.ipynb
index 79b779a96a..c3f9ec0fcb 100644
--- a/docs/notebooks/11-extract-adjoint-solutions.ipynb
+++ b/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",
diff --git a/docs/source/interpolation.rst b/docs/source/interpolation.rst
index b079a9e352..66082da969 100644
--- a/docs/source/interpolation.rst
+++ b/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
diff --git a/firedrake/__future__.py b/firedrake/__future__.py
index b0c67effdb..d3656f79ad 100644
--- a/firedrake/__future__.py
+++ b/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)
diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py
index 351b393394..ccf0a81cfc 100644
--- a/firedrake/interpolation.py
+++ b/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
 
diff --git a/firedrake/preconditioners/hypre_ads.py b/firedrake/preconditioners/hypre_ads.py
index 13dc0d1a50..b639c21f05 100644
--- a/firedrake/preconditioners/hypre_ads.py
+++ b/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
 
diff --git a/firedrake/preconditioners/hypre_ams.py b/firedrake/preconditioners/hypre_ams.py
index 9616b142b9..5f995c3d90 100644
--- a/firedrake/preconditioners/hypre_ams.py
+++ b/firedrake/preconditioners/hypre_ams.py
@@ -1,4 +1,4 @@
-import firedrake.assemble as assemble
+import firedrake as fd
 from firedrake.preconditioners.base import PCBase
 from firedrake.petsc import PETSc
 from firedrake.functionspace import FunctionSpace, VectorFunctionSpace
@@ -6,9 +6,6 @@
 from firedrake.dmhooks import get_function_space
 from firedrake.utils import complex_mode
 from firedrake_citations import Citations
-from firedrake import SpatialCoordinate
-from firedrake.__future__ import interpolate
-from ufl import grad
 from pyop2.utils import as_tuple
 
 __all__ = ("HypreAMS",)
@@ -50,7 +47,7 @@ def initialize(self, obj):
         P1 = FunctionSpace(mesh, "Lagrange", 1)
         G_callback = appctx.get("get_gradient", None)
         if G_callback is None:
-            G = chop(assemble.assemble(interpolate(grad(TestFunction(P1)), V)).petscmat)
+            G = chop(fd.assemble(fd.interpolate(fd.grad(TestFunction(P1)), V)).petscmat)
         else:
             G = G_callback(P1, V)
 
@@ -68,8 +65,8 @@ def initialize(self, obj):
             pc.setHYPRESetBetaPoissonMatrix(None)
 
         VectorP1 = VectorFunctionSpace(mesh, "Lagrange", 1)
-        interp = interpolate(SpatialCoordinate(mesh), VectorP1)
-        pc.setCoordinates(assemble.assemble(interp).dat.data_ro.copy())
+        interp = fd.interpolate(fd.SpatialCoordinate(mesh), VectorP1)
+        pc.setCoordinates(fd.assemble(interp).dat.data_ro.copy())
         pc.setFromOptions()
         self.pc = pc
 
diff --git a/tests/firedrake/adjoint/test_assignment.py b/tests/firedrake/adjoint/test_assignment.py
index c4b9fcd467..93bee8f265 100644
--- a/tests/firedrake/adjoint/test_assignment.py
+++ b/tests/firedrake/adjoint/test_assignment.py
@@ -1,7 +1,6 @@
 import pytest
 
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.adjoint import *
 
 from numpy.random import rand
diff --git a/tests/firedrake/adjoint/test_disk_checkpointing.py b/tests/firedrake/adjoint/test_disk_checkpointing.py
index 6c4e5deb45..f264e5cb11 100644
--- a/tests/firedrake/adjoint/test_disk_checkpointing.py
+++ b/tests/firedrake/adjoint/test_disk_checkpointing.py
@@ -1,7 +1,6 @@
 import pytest
 
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.adjoint import *
 from firedrake.adjoint_utils.checkpointing import disk_checkpointing
 from checkpoint_schedules import SingleDiskStorageSchedule
diff --git a/tests/firedrake/adjoint/test_projection.py b/tests/firedrake/adjoint/test_projection.py
index 1b38b1170f..d77e6a9df0 100644
--- a/tests/firedrake/adjoint/test_projection.py
+++ b/tests/firedrake/adjoint/test_projection.py
@@ -1,7 +1,6 @@
 import pytest
 
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.adjoint import *
 
 from numpy.random import rand
diff --git a/tests/firedrake/adjoint/test_shape_derivatives.py b/tests/firedrake/adjoint/test_shape_derivatives.py
index 88cee34f78..cb23522ce5 100644
--- a/tests/firedrake/adjoint/test_shape_derivatives.py
+++ b/tests/firedrake/adjoint/test_shape_derivatives.py
@@ -2,7 +2,6 @@
 
 import numpy as np
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.adjoint import *
 from pyadjoint import taylor_to_dict
 
diff --git a/tests/firedrake/extrusion/test_variable_layers_bcs.py b/tests/firedrake/extrusion/test_variable_layers_bcs.py
index 6c7ceb4b55..1dc3498415 100644
--- a/tests/firedrake/extrusion/test_variable_layers_bcs.py
+++ b/tests/firedrake/extrusion/test_variable_layers_bcs.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.utils import IntType
 import pytest
 import numpy
diff --git a/tests/firedrake/extrusion/test_variable_layers_numbering.py b/tests/firedrake/extrusion/test_variable_layers_numbering.py
index 326a06a499..e13358f57d 100644
--- a/tests/firedrake/extrusion/test_variable_layers_numbering.py
+++ b/tests/firedrake/extrusion/test_variable_layers_numbering.py
@@ -1,7 +1,6 @@
 import pytest
 import numpy
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.mesh import plex_from_cell_list
 from firedrake.utils import IntType
 
diff --git a/tests/firedrake/extrusion/test_variable_layers_poisson.py b/tests/firedrake/extrusion/test_variable_layers_poisson.py
index 7c3244eb53..7ab9bdac21 100644
--- a/tests/firedrake/extrusion/test_variable_layers_poisson.py
+++ b/tests/firedrake/extrusion/test_variable_layers_poisson.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 import numpy
 from firedrake.utils import IntType
 
diff --git a/tests/firedrake/extrusion/test_variable_layers_steady_advection.py b/tests/firedrake/extrusion/test_variable_layers_steady_advection.py
index f26838bfb9..d8439b73bd 100644
--- a/tests/firedrake/extrusion/test_variable_layers_steady_advection.py
+++ b/tests/firedrake/extrusion/test_variable_layers_steady_advection.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 import numpy
 from firedrake.utils import IntType
 
diff --git a/tests/firedrake/macro/test_macro_multigrid.py b/tests/firedrake/macro/test_macro_multigrid.py
index 8f880bd48d..e3d8b11e9f 100644
--- a/tests/firedrake/macro/test_macro_multigrid.py
+++ b/tests/firedrake/macro/test_macro_multigrid.py
@@ -1,7 +1,6 @@
 import numpy
 import pytest
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.utils import complex_mode
 
 
diff --git a/tests/firedrake/multigrid/test_grid_transfer.py b/tests/firedrake/multigrid/test_grid_transfer.py
index 40b21bb3b2..09c279baec 100644
--- a/tests/firedrake/multigrid/test_grid_transfer.py
+++ b/tests/firedrake/multigrid/test_grid_transfer.py
@@ -1,7 +1,6 @@
 import pytest
 import numpy
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.utils import complex_mode
 
 
diff --git a/tests/firedrake/output/test_plotting.py b/tests/firedrake/output/test_plotting.py
index b3195d0038..2cf840e2e9 100644
--- a/tests/firedrake/output/test_plotting.py
+++ b/tests/firedrake/output/test_plotting.py
@@ -1,7 +1,6 @@
 import pytest
 import numpy as np
 from firedrake import *
-from firedrake.__future__ import *
 
 try:
     from firedrake.pyplot import *
diff --git a/tests/firedrake/regression/test_adjoint_operators.py b/tests/firedrake/regression/test_adjoint_operators.py
index b7ecdd6ed0..6410d74957 100644
--- a/tests/firedrake/regression/test_adjoint_operators.py
+++ b/tests/firedrake/regression/test_adjoint_operators.py
@@ -6,7 +6,6 @@
 
 from firedrake import *
 from firedrake.adjoint import *
-from firedrake.__future__ import *
 
 
 @pytest.fixture(autouse=True)
diff --git a/tests/firedrake/regression/test_cellorigin.py b/tests/firedrake/regression/test_cellorigin.py
index c763de56b0..754fbf23d3 100644
--- a/tests/firedrake/regression/test_cellorigin.py
+++ b/tests/firedrake/regression/test_cellorigin.py
@@ -2,7 +2,6 @@
 import numpy as np
 from ufl.geometry import CellOrigin
 from firedrake import *
-from firedrake.__future__ import *
 
 
 @pytest.fixture(params=["interval", "triangle", "quadrilateral", "tetrahedron"])
diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py
index 04d9254963..ae3993534c 100644
--- a/tests/firedrake/regression/test_interp_dual.py
+++ b/tests/firedrake/regression/test_interp_dual.py
@@ -1,7 +1,6 @@
 import pytest
 import numpy as np
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.utils import complex_mode
 import ufl
 
diff --git a/tests/firedrake/regression/test_interpolate.py b/tests/firedrake/regression/test_interpolate.py
index e8860a32e7..45a2d9c89a 100644
--- a/tests/firedrake/regression/test_interpolate.py
+++ b/tests/firedrake/regression/test_interpolate.py
@@ -2,7 +2,6 @@
 import numpy as np
 import pytest
 from firedrake import *
-from firedrake.__future__ import *
 
 cwd = abspath(dirname(__file__))
 
diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py
index c19b028a19..c3e3c39ece 100644
--- a/tests/firedrake/regression/test_interpolate_cross_mesh.py
+++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.petsc import DEFAULT_PARTITIONER
 from firedrake.ufl_expr import extract_unique_domain
 import numpy as np
diff --git a/tests/firedrake/regression/test_interpolate_p3intmoments.py b/tests/firedrake/regression/test_interpolate_p3intmoments.py
index 6ef0ef34d2..5a7e18c1d6 100644
--- a/tests/firedrake/regression/test_interpolate_p3intmoments.py
+++ b/tests/firedrake/regression/test_interpolate_p3intmoments.py
@@ -1,7 +1,6 @@
 import pytest
 import numpy as np
 from firedrake import *
-from firedrake.__future__ import *
 from FIAT.dual_set import DualSet
 from FIAT.finite_element import CiarletElement
 from FIAT.reference_element import UFCInterval
diff --git a/tests/firedrake/regression/test_interpolate_vs_project.py b/tests/firedrake/regression/test_interpolate_vs_project.py
index 84fabd6a45..3a0fcc571f 100644
--- a/tests/firedrake/regression/test_interpolate_vs_project.py
+++ b/tests/firedrake/regression/test_interpolate_vs_project.py
@@ -1,7 +1,6 @@
 import numpy as np
 import pytest
 from firedrake import *
-from firedrake.__future__ import *
 
 
 @pytest.fixture(params=["square", "cube"], scope="module")
diff --git a/tests/firedrake/regression/test_interpolation_manual.py b/tests/firedrake/regression/test_interpolation_manual.py
index ffe515ddd1..3b2dcfa041 100644
--- a/tests/firedrake/regression/test_interpolation_manual.py
+++ b/tests/firedrake/regression/test_interpolation_manual.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 import pytest
 import numpy as np
 
diff --git a/tests/firedrake/regression/test_interpolation_nodes.py b/tests/firedrake/regression/test_interpolation_nodes.py
index 6c6e338afb..13b52aa193 100644
--- a/tests/firedrake/regression/test_interpolation_nodes.py
+++ b/tests/firedrake/regression/test_interpolation_nodes.py
@@ -1,7 +1,6 @@
 import numpy as np
 import pytest
 from firedrake import *
-from firedrake.__future__ import *
 
 '''
 The spaces N1div, N1curl, N2div and N2curl have the special property that the interpolation in these
diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py
index 6e728ef653..51371d36da 100644
--- a/tests/firedrake/regression/test_netgen.py
+++ b/tests/firedrake/regression/test_netgen.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import interpolate
 import numpy as np
 import pytest
 
diff --git a/tests/firedrake/regression/test_norm.py b/tests/firedrake/regression/test_norm.py
index d66c55ee85..35762c4c38 100644
--- a/tests/firedrake/regression/test_norm.py
+++ b/tests/firedrake/regression/test_norm.py
@@ -1,6 +1,5 @@
 import pytest
 from firedrake import *
-from firedrake.__future__ import *
 import numpy
 
 
diff --git a/tests/firedrake/regression/test_nullspace.py b/tests/firedrake/regression/test_nullspace.py
index e4f436ff8d..5b3538c41e 100644
--- a/tests/firedrake/regression/test_nullspace.py
+++ b/tests/firedrake/regression/test_nullspace.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.petsc import PETSc
 import pytest
 import numpy as np
diff --git a/tests/firedrake/regression/test_octahedral_hemisphere.py b/tests/firedrake/regression/test_octahedral_hemisphere.py
index be2c0696c9..da8ed1aada 100644
--- a/tests/firedrake/regression/test_octahedral_hemisphere.py
+++ b/tests/firedrake/regression/test_octahedral_hemisphere.py
@@ -1,7 +1,6 @@
 from __future__ import absolute_import, print_function, division
 import pytest
 from firedrake import *
-from firedrake.__future__ import *
 import numpy
 
 
diff --git a/tests/firedrake/regression/test_par_loops.py b/tests/firedrake/regression/test_par_loops.py
index 8222e0d5c7..f4df8f51fa 100644
--- a/tests/firedrake/regression/test_par_loops.py
+++ b/tests/firedrake/regression/test_par_loops.py
@@ -2,7 +2,6 @@
 import numpy as np
 from firedrake import *
 from firedrake.utils import ScalarType_c
-from firedrake.__future__ import *
 
 
 @pytest.fixture(scope="module")
diff --git a/tests/firedrake/regression/test_project_interp_KMV.py b/tests/firedrake/regression/test_project_interp_KMV.py
index 44ca5ec59c..7eef071841 100644
--- a/tests/firedrake/regression/test_project_interp_KMV.py
+++ b/tests/firedrake/regression/test_project_interp_KMV.py
@@ -1,7 +1,6 @@
 import pytest
 import numpy as np
 from firedrake import *
-from firedrake.__future__ import *
 import finat
 
 
diff --git a/tests/firedrake/regression/test_real_space.py b/tests/firedrake/regression/test_real_space.py
index b989ae2c86..74f4015c66 100644
--- a/tests/firedrake/regression/test_real_space.py
+++ b/tests/firedrake/regression/test_real_space.py
@@ -2,7 +2,6 @@
 import numpy as np
 
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.petsc import DEFAULT_DIRECT_SOLVER
 
 
diff --git a/tests/firedrake/slate/test_facet_tensors.py b/tests/firedrake/slate/test_facet_tensors.py
index 834db4e09b..eacf020266 100644
--- a/tests/firedrake/slate/test_facet_tensors.py
+++ b/tests/firedrake/slate/test_facet_tensors.py
@@ -1,7 +1,6 @@
 import pytest
 import numpy as np
 from firedrake import *
-from firedrake.__future__ import *
 
 
 @pytest.fixture(scope='module', params=[False, True])
diff --git a/tests/firedrake/slate/test_facet_tensors_extr.py b/tests/firedrake/slate/test_facet_tensors_extr.py
index fb4443ac76..1a0d058b67 100644
--- a/tests/firedrake/slate/test_facet_tensors_extr.py
+++ b/tests/firedrake/slate/test_facet_tensors_extr.py
@@ -1,7 +1,6 @@
 import pytest
 import numpy as np
 from firedrake import *
-from firedrake.__future__ import *
 
 
 @pytest.fixture(scope='module', params=[False, True])
diff --git a/tests/firedrake/vertexonly/test_interpolation_from_parent.py b/tests/firedrake/vertexonly/test_interpolation_from_parent.py
index c722d7b9b3..275c729cf2 100644
--- a/tests/firedrake/vertexonly/test_interpolation_from_parent.py
+++ b/tests/firedrake/vertexonly/test_interpolation_from_parent.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 import pytest
 import numpy as np
 from functools import reduce
diff --git a/tests/firedrake/vertexonly/test_point_eval_immersed_manifold.py b/tests/firedrake/vertexonly/test_point_eval_immersed_manifold.py
index 30d8057748..bfae5673e1 100644
--- a/tests/firedrake/vertexonly/test_point_eval_immersed_manifold.py
+++ b/tests/firedrake/vertexonly/test_point_eval_immersed_manifold.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 import numpy as np
 import pytest
 
diff --git a/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py b/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py
index df28b196f2..9975aa6dd1 100644
--- a/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py
+++ b/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py
@@ -1,7 +1,6 @@
 import pytest
 import numpy as np
 from firedrake import *
-from firedrake.__future__ import *
 from pyadjoint.tape import get_working_tape, pause_annotation
 
 
diff --git a/tests/firedrake/vertexonly/test_swarm.py b/tests/firedrake/vertexonly/test_swarm.py
index 4bd5cde939..555351cd01 100644
--- a/tests/firedrake/vertexonly/test_swarm.py
+++ b/tests/firedrake/vertexonly/test_swarm.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 from firedrake.utils import IntType, RealType
 import pytest
 import numpy as np
diff --git a/tests/firedrake/vertexonly/test_vertex_only_fs.py b/tests/firedrake/vertexonly/test_vertex_only_fs.py
index ce071267b6..6fa9fddd7c 100644
--- a/tests/firedrake/vertexonly/test_vertex_only_fs.py
+++ b/tests/firedrake/vertexonly/test_vertex_only_fs.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 import pytest
 import numpy as np
 from mpi4py import MPI
diff --git a/tests/firedrake/vertexonly/test_vertex_only_manual.py b/tests/firedrake/vertexonly/test_vertex_only_manual.py
index 7578b5c69b..f44b30758e 100644
--- a/tests/firedrake/vertexonly/test_vertex_only_manual.py
+++ b/tests/firedrake/vertexonly/test_vertex_only_manual.py
@@ -1,5 +1,4 @@
 from firedrake import *
-from firedrake.__future__ import *
 import pytest