refactor interpolation

demos/boussinesq/boussinesq.py.rst

@@ -184,7 +184,7 @@ implements a boundary condition that fixes a field at a single point. ::
 
            # Take the basis function with the largest abs value at bc_point
            v = TestFunction(V)
-           F = assemble(Interpolate(inner(v, v), Fvom))
+           F = assemble(interpolate(inner(v, v), Fvom))
            with F.dat.vec as Fvec:
                max_index, _ = Fvec.max()
            nodes = V.dof_dset.lgmap.applyInverse([max_index])

demos/multicomponent/multicomponent.py.rst

@@ -521,7 +521,7 @@ mathematically valid to do this)::
 
            # Take the basis function with the largest abs value at bc_point
            v = TestFunction(V)
-           F = assemble(Interpolate(inner(v, v), Fvom))
+           F = assemble(interpolate(inner(v, v), Fvom))
            with F.dat.vec as Fvec:
                max_index, _ = Fvec.max()
            nodes = V.dof_dset.lgmap.applyInverse([max_index])

firedrake/assemble.py

@@ -23,6 +23,7 @@
 from firedrake.bcs import DirichletBC, EquationBC, EquationBCSplit
 from firedrake.functionspaceimpl import WithGeometry, FunctionSpace, FiredrakeDualSpace
 from firedrake.functionspacedata import entity_dofs_key, entity_permutations_key
+from firedrake.interpolation import get_interpolator
 from firedrake.petsc import PETSc
 from firedrake.slate import slac, slate
 from firedrake.slate.slac.kernel_builder import CellFacetKernelArg, LayerCountKernelArg
@@ -578,17 +579,9 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
             rank = len(expr.arguments())
             if rank > 2:
                 raise ValueError("Cannot assemble an Interpolate with more than two arguments")
-            # Get the target space
-            V = v.function_space().dual()
-
-            # Get the interpolator
-            interp_data = expr.interp_data.copy()
-            default_missing_val = interp_data.pop('default_missing_val', None)
-            if rank == 1 and isinstance(tensor, firedrake.Function):
-                V = tensor
-            interpolator = firedrake.Interpolator(expr, V, **interp_data)
-            # Assembly
-            return interpolator.assemble(tensor=tensor, default_missing_val=default_missing_val)
+
+            interpolator = get_interpolator(expr)
+            return interpolator.assemble(tensor=tensor)
         elif tensor and isinstance(expr, (firedrake.Function, firedrake.Cofunction, firedrake.MatrixBase)):
             return tensor.assign(expr)
         elif tensor and isinstance(expr, ufl.ZeroBaseForm):

firedrake/bcs.py

@@ -1,7 +1,7 @@
 # A module implementing strong (Dirichlet) boundary conditions.
 import numpy as np
 
-import functools
+from functools import partial, reduce
 import itertools
 
 import ufl
@@ -167,7 +167,7 @@ def hermite_stride(bcnodes):
                     # Edge conditions have only been tested with Lagrange elements.
                     # Need to expand the list.
                     bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
-                bcnodes1 = functools.reduce(np.intersect1d, bcnodes1)
+                bcnodes1 = reduce(np.intersect1d, bcnodes1)
                 bcnodes.append(bcnodes1)
         return np.concatenate(bcnodes)
 
@@ -339,6 +339,7 @@ def reconstruct(self, field=None, V=None, g=None, sub_domain=None, use_split=Fal
     @function_arg.setter
     def function_arg(self, g):
         '''Set the value of this boundary condition.'''
+        from firedrake.interpolation import get_interpolator
         try:
             # Clear any previously set update function
             del self._function_arg_update
@@ -359,12 +360,11 @@ def function_arg(self, g):
                 raise RuntimeError(f"Provided boundary value {g} does not match shape of space")
             try:
                 self._function_arg = firedrake.Function(V)
-                # Use `Interpolator` instead of assembling an `Interpolate` form
-                # as the expression compilation needs to happen at this stage to
-                # determine if we should use interpolation or projection
-                #  -> e.g. interpolation may not be supported for the element.
-                self._function_arg_update = firedrake.Interpolator(g, self._function_arg)._interpolate
-            except (NotImplementedError, AttributeError):
+                interpolator = get_interpolator(firedrake.interpolate(g, V))
+                # Call this here to check if the element supports interpolation
+                interpolator._build_callable()
+                self._function_arg_update = partial(interpolator.assemble, tensor=self._function_arg)
+            except (ValueError, NotImplementedError):
                 # Element doesn't implement interpolation
                 self._function_arg = firedrake.Function(V).project(g)
                 self._function_arg_update = firedrake.Projector(g, self._function_arg).project

firedrake/cofunction.py

@@ -318,7 +318,7 @@ def interpolate(self,
         Parameters
         ----------
         expression
-            A dual UFL expression to interpolate.
+            A UFL BaseForm to adjoint interpolate.
         ad_block_tag
             An optional string for tagging the resulting assemble
             block on the Pyadjoint tape.
@@ -331,9 +331,9 @@ def interpolate(self,
         firedrake.cofunction.Cofunction
             Returns `self`
         """
-        from firedrake import interpolation, assemble
+        from firedrake import interpolate, assemble
         v, = self.arguments()
-        interp = interpolation.Interpolate(v, expression, **kwargs)
+        interp = interpolate(v, expression, **kwargs)
         return assemble(interp, tensor=self, ad_block_tag=ad_block_tag)
 
     @property

firedrake/external_operators/point_expr_operator.py

@@ -5,7 +5,7 @@
 
 import firedrake.ufl_expr as ufl_expr
 from firedrake.assemble import assemble
-from firedrake.interpolation import Interpolate
+from firedrake.interpolation import interpolate
 from firedrake.external_operators import AbstractExternalOperator, assemble_method
 
 
@@ -58,7 +58,7 @@ def assemble_operator(self, *args, **kwargs):
         V = self.function_space()
         expr = as_ufl(self.expr(*self.ufl_operands))
         if len(V) < 2:
-            interp = Interpolate(expr, self.function_space())
+            interp = interpolate(expr, self.function_space())
             return assemble(interp)
         # Interpolation of UFL expressions for mixed functions is not yet supported
         # -> `Function.assign` might be enough in some cases.
@@ -72,7 +72,7 @@ def assemble_operator(self, *args, **kwargs):
     def assemble_Jacobian_action(self, *args, **kwargs):
         V = self.function_space()
         expr = as_ufl(self.expr(*self.ufl_operands))
-        interp = Interpolate(expr, V)
+        interp = interpolate(expr, V)
 
         u, = [e for i, e in enumerate(self.ufl_operands) if self.derivatives[i] == 1]
         w = self.argument_slots()[-1]
@@ -83,7 +83,7 @@ def assemble_Jacobian_action(self, *args, **kwargs):
     def assemble_Jacobian(self, *args, assembly_opts, **kwargs):
         V = self.function_space()
         expr = as_ufl(self.expr(*self.ufl_operands))
-        interp = Interpolate(expr, V)
+        interp = interpolate(expr, V)
 
         u, = [e for i, e in enumerate(self.ufl_operands) if self.derivatives[i] == 1]
         jac = ufl_expr.derivative(interp, u)
@@ -99,7 +99,7 @@ def assemble_Jacobian_adjoint(self, *args, assembly_opts, **kwargs):
     def assemble_Jacobian_adjoint_action(self, *args, **kwargs):
         V = self.function_space()
         expr = as_ufl(self.expr(*self.ufl_operands))
-        interp = Interpolate(expr, V)
+        interp = interpolate(expr, V)
 
         u, = [e for i, e in enumerate(self.ufl_operands) if self.derivatives[i] == 1]
         ustar = self.argument_slots()[0]

firedrake/function.py

@@ -382,9 +382,9 @@ def interpolate(self,
         firedrake.function.Function
             Returns `self`
         """
-        from firedrake import interpolation, assemble
+        from firedrake import interpolate, assemble
         V = self.function_space()
-        interp = interpolation.Interpolate(expression, V, **kwargs)
+        interp = interpolate(expression, V, **kwargs)
         return assemble(interp, tensor=self, ad_block_tag=ad_block_tag)
 
     def zero(self, subset=None):
@@ -701,7 +701,7 @@ def __init__(self, domain, point):
         self.point = point
 
     def __str__(self):
-        return "domain %s does not contain point %s" % (self.domain, self.point)
+        return f"Domain {self.domain} does not contain point {self.point}"
 
 
 class PointEvaluator:
@@ -716,7 +716,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             The mesh on which to embed the points.
         points : numpy.ndarray | list
             Array or list of points to evaluate at.
-        tolerance : float | None
+        tolerance : Optional[float]
             Tolerance to use when checking if a point is in a cell.
             If ``None`` (the default), the ``tolerance`` of the ``mesh`` is used.
         missing_points_behaviour : str