From d286b36a6475cf147dd17aa533e149c464a969d1 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 23 Sep 2025 15:51:07 +0100 Subject: [PATCH 01/21] simplify . --- firedrake/interpolation.py | 30 ++++-------------------------- 1 file changed, 4 insertions(+), 26 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 4fb288fb4d..04b0993f23 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -174,32 +174,10 @@ def interpolate(expr, V, subset=None, access=None, allow_missing_dofs=False, def reduction (hence using MIN will compute the MIN between the existing values and any new values). """ - if isinstance(V, (Cofunction, Coargument)): - dual_arg = V - elif isinstance(V, ufl.BaseForm): - rank = len(V.arguments()) - if rank == 1: - dual_arg = V - else: - raise TypeError(f"Expected a one-form, provided form had {rank} arguments") - elif isinstance(V, functionspaceimpl.WithGeometry): - dual_arg = Coargument(V.dual(), 0) - expr_args = extract_arguments(ufl.as_ufl(expr)) - if expr_args and expr_args[0].number() == 0: - warnings.warn("Passing argument numbered 0 in expression for forward interpolation is deprecated. " - "Use a TrialFunction in the expression.") - v, = expr_args - expr = replace(expr, {v: v.reconstruct(number=1)}) - else: - raise TypeError(f"V must be a FunctionSpace, Cofunction, Coargument or one-form, not a {type(V).__name__}") - - interp = Interpolate(expr, dual_arg, - subset=subset, access=access, - allow_missing_dofs=allow_missing_dofs, - default_missing_val=default_missing_val, - matfree=matfree) - - return interp + return Interpolate( + expr, V, subset=subset, access=access, allow_missing_dofs=allow_missing_dofs, + default_missing_val=default_missing_val, matfree=matfree + ) class Interpolator(abc.ABC): From d930a1b2697472d755d88267197b95197954ef70 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 23 Sep 2025 15:55:12 +0100 Subject: [PATCH 02/21] tidy function interpolate --- firedrake/function.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/firedrake/function.py b/firedrake/function.py index a628ac6599..fed408c8cc 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -384,7 +384,10 @@ def interpolate(self, """ from firedrake import interpolation, assemble V = self.function_space() - interp = interpolation.Interpolate(expression, V, **kwargs) + interp = interpolate( + expression, V, subset=subset, allow_missing_dofs=allow_missing_dofs, + default_missing_val=default_missing_val + ) return assemble(interp, tensor=self, ad_block_tag=ad_block_tag) def zero(self, subset=None): From 7f510136d443826ee9fb241b2e1e6b8353cdc9b6 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 23 Sep 2025 16:12:40 +0100 Subject: [PATCH 03/21] create Coargument in Firedrake --- firedrake/interpolation.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 04b0993f23..deb0df1f50 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -174,6 +174,11 @@ def interpolate(expr, V, subset=None, access=None, allow_missing_dofs=False, def reduction (hence using MIN will compute the MIN between the existing values and any new values). """ + if isinstance(V, functionspaceimpl.WithGeometry): + # Need to create a Firedrake Argument so that it has a .function_space() method + expr_args = extract_arguments(ufl.as_ufl(expr)) + is_adjoint = len(expr_args) and expr_args[0].number() == 0 + V = Argument(V.dual(), 1 if is_adjoint else 0) return Interpolate( expr, V, subset=subset, access=access, allow_missing_dofs=allow_missing_dofs, default_missing_val=default_missing_val, matfree=matfree From c2f5476abfea317ec99b6c30980ac39016bdf60b Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 23 Sep 2025 17:28:25 +0100 Subject: [PATCH 04/21] Change `Interpolate` to `interpolate` --- demos/boussinesq/boussinesq.py.rst | 2 +- demos/multicomponent/multicomponent.py.rst | 2 +- firedrake/external_operators/point_expr_operator.py | 10 +++++----- firedrake/interpolation.py | 6 +++--- firedrake/mesh.py | 2 +- firedrake/mg/utils.py | 2 +- firedrake/preconditioners/gtmg.py | 4 ++-- firedrake/preconditioners/patch.py | 4 ++-- firedrake/pyplot/mpl.py | 12 ++++++------ firedrake/utility_meshes.py | 8 ++++---- tests/firedrake/adjoint/test_reduced_functional.py | 4 ++-- .../external_operators/test_external_operators.py | 2 +- tests/firedrake/multigrid/test_poisson_gtmg.py | 2 +- tests/firedrake/regression/test_adjoint_operators.py | 2 +- tests/firedrake/regression/test_interpolate_zany.py | 6 +++--- tests/firedrake/submesh/test_submesh_interpolate.py | 4 ++-- 16 files changed, 36 insertions(+), 36 deletions(-) diff --git a/demos/boussinesq/boussinesq.py.rst b/demos/boussinesq/boussinesq.py.rst index edfdf3c1a5..5cc6708cf0 100644 --- a/demos/boussinesq/boussinesq.py.rst +++ b/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]) diff --git a/demos/multicomponent/multicomponent.py.rst b/demos/multicomponent/multicomponent.py.rst index bf74e5d2e0..7a29d6ef1c 100644 --- a/demos/multicomponent/multicomponent.py.rst +++ b/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]) diff --git a/firedrake/external_operators/point_expr_operator.py b/firedrake/external_operators/point_expr_operator.py index 3aa40e1d5b..4e7183e47f 100644 --- a/firedrake/external_operators/point_expr_operator.py +++ b/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] diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index deb0df1f50..bf3fcb5013 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -511,7 +511,7 @@ def __init__( from firedrake.assemble import assemble V_dest_vec = firedrake.VectorFunctionSpace(dest_mesh, ufl_scalar_element) - f_dest_node_coords = Interpolate(dest_mesh.coordinates, V_dest_vec) + f_dest_node_coords = interpolate(dest_mesh.coordinates, V_dest_vec) f_dest_node_coords = assemble(f_dest_node_coords) dest_node_coords = f_dest_node_coords.dat.data_ro.reshape(-1, dest_mesh_gdim) try: @@ -536,7 +536,7 @@ def __init__( else: fs_type = partial(firedrake.TensorFunctionSpace, shape=shape) P0DG_vom = fs_type(self.vom_dest_node_coords_in_src_mesh, "DG", 0) - self.point_eval_interpolate = Interpolate(self.expr_renumbered, P0DG_vom) + self.point_eval_interpolate = interpolate(self.expr_renumbered, P0DG_vom) # The parallel decomposition of the nodes of V_dest in the DESTINATION # mesh (dest_mesh) is retrieved using the input_ordering attribute of the # VOM. This again is an interpolation operation, which, under the hood @@ -544,7 +544,7 @@ def __init__( P0DG_vom_i_o = fs_type( self.vom_dest_node_coords_in_src_mesh.input_ordering, "DG", 0 ) - self.to_input_ordering_interpolate = Interpolate( + self.to_input_ordering_interpolate = interpolate( firedrake.TrialFunction(P0DG_vom), P0DG_vom_i_o ) # The P0DG function outputted by the above interpolation has the diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 87757a2e54..14ebf94d0e 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -4140,7 +4140,7 @@ def _parent_mesh_embedding( # nessesary, to other processes. P0DG = functionspace.FunctionSpace(parent_mesh, "DG", 0) with stop_annotating(): - visible_ranks = interpolation.Interpolate( + visible_ranks = interpolation.interpolate( constant.Constant(parent_mesh.comm.rank), P0DG ) visible_ranks = assemble(visible_ranks).dat.data_ro_with_halos.real diff --git a/firedrake/mg/utils.py b/firedrake/mg/utils.py index 37832b64dc..886cc7530c 100644 --- a/firedrake/mg/utils.py +++ b/firedrake/mg/utils.py @@ -143,7 +143,7 @@ def physical_node_locations(V): Vc = V.collapse().reconstruct(element=finat.ufl.VectorElement(element, dim=mesh.geometric_dimension())) # FIXME: This is unsafe for DG coordinates and CG target spaces. - locations = firedrake.assemble(firedrake.Interpolate(firedrake.SpatialCoordinate(mesh), Vc)) + locations = firedrake.assemble(firedrake.interpolate(firedrake.SpatialCoordinate(mesh), Vc)) return cache.setdefault(key, locations) diff --git a/firedrake/preconditioners/gtmg.py b/firedrake/preconditioners/gtmg.py index 6ce73cd6b4..2ac5df9a5d 100644 --- a/firedrake/preconditioners/gtmg.py +++ b/firedrake/preconditioners/gtmg.py @@ -4,7 +4,7 @@ from firedrake.petsc import PETSc from firedrake.preconditioners.base import PCBase from firedrake.parameters import parameters -from firedrake.interpolation import Interpolate +from firedrake.interpolation import interpolate from firedrake.solving_utils import _SNESContext from firedrake.matrix_free.operators import ImplicitMatrixContext import firedrake.dmhooks as dmhooks @@ -155,7 +155,7 @@ def initialize(self, pc): # Create interpolation matrix from coarse space to fine space fine_space = ctx.J.arguments()[0].function_space() coarse_test, coarse_trial = coarse_operator.arguments() - interp = assemble(Interpolate(coarse_trial, fine_space)) + interp = assemble(interpolate(coarse_trial, fine_space)) interp_petscmat = interp.petscmat restr_petscmat = appctx.get("restriction_matrix", None) diff --git a/firedrake/preconditioners/patch.py b/firedrake/preconditioners/patch.py index c139c07f7a..2fbe469e0f 100644 --- a/firedrake/preconditioners/patch.py +++ b/firedrake/preconditioners/patch.py @@ -4,7 +4,7 @@ from firedrake.solving_utils import _SNESContext from firedrake.utils import cached_property, complex_mode, IntType from firedrake.dmhooks import get_appctx, push_appctx, pop_appctx -from firedrake.interpolation import Interpolate +from firedrake.interpolation import interpolate from collections import namedtuple import operator @@ -660,7 +660,7 @@ def sort_entities(self, dm, axis, dir, ndiv=None, divisions=None): # with access descriptor MAX to define a consistent opinion # about where the vertices are. CGk = V.reconstruct(family="Lagrange") - coordinates = assemble(Interpolate(coordinates, CGk, access=op2.MAX)) + coordinates = assemble(interpolate(coordinates, CGk, access=op2.MAX)) select = partial(select_entity, dm=dm, exclude="pyop2_ghost") entities = [(p, self.coords(dm, p, coordinates)) for p in diff --git a/firedrake/pyplot/mpl.py b/firedrake/pyplot/mpl.py index 3cf010a1c9..d6a7aa5112 100644 --- a/firedrake/pyplot/mpl.py +++ b/firedrake/pyplot/mpl.py @@ -18,7 +18,7 @@ import mpl_toolkits.mplot3d from mpl_toolkits.mplot3d.art3d import Line3DCollection, Poly3DCollection from math import factorial -from firedrake import (Interpolate, sqrt, inner, Function, SpatialCoordinate, +from firedrake import (interpolate, sqrt, inner, Function, SpatialCoordinate, FunctionSpace, VectorFunctionSpace, PointNotInDomainError, Constant, assemble, dx) from firedrake.mesh import MeshGeometry @@ -120,7 +120,7 @@ def triplot(mesh, axes=None, interior_kw={}, boundary_kw={}): if element.degree() != 1: # Interpolate to piecewise linear. V = VectorFunctionSpace(mesh, element.family(), 1) - coordinates = assemble(Interpolate(coordinates, V)) + coordinates = assemble(interpolate(coordinates, V)) coords = toreal(coordinates.dat.data_ro_with_halos, "real") result = [] @@ -215,7 +215,7 @@ def _plot_2d_field(method_name, function, *args, complex_component="real", **kwa if len(function.ufl_shape) == 1: element = function.ufl_element().sub_elements[0] Q = FunctionSpace(mesh, element) - function = assemble(Interpolate(sqrt(inner(function, function)), Q)) + function = assemble(interpolate(sqrt(inner(function, function)), Q)) num_sample_points = kwargs.pop("num_sample_points", 10) function_plotter = FunctionPlotter(mesh, num_sample_points) @@ -326,7 +326,7 @@ def trisurf(function, *args, complex_component="real", **kwargs): if len(function.ufl_shape) == 1: element = function.ufl_element().sub_elements[0] Q = FunctionSpace(mesh, element) - function = assemble(Interpolate(sqrt(inner(function, function)), Q)) + function = assemble(interpolate(sqrt(inner(function, function)), Q)) num_sample_points = kwargs.pop("num_sample_points", 10) function_plotter = FunctionPlotter(mesh, num_sample_points) @@ -355,7 +355,7 @@ def quiver(function, *, complex_component="real", **kwargs): coords = toreal(extract_unique_domain(function).coordinates.dat.data_ro, "real") V = extract_unique_domain(function).coordinates.function_space() - function_interp = assemble(Interpolate(function, V)) + function_interp = assemble(interpolate(function, V)) vals = toreal(function_interp.dat.data_ro, complex_component) C = np.linalg.norm(vals, axis=1) return axes.quiver(*(coords.T), *(vals.T), C, **kwargs) @@ -816,7 +816,7 @@ def _bezier_plot(function, axes, complex_component="real", **kwargs): mesh = function.function_space().mesh() if deg == 0: V = FunctionSpace(mesh, "DG", 1) - interp = assemble(Interpolate(function, V)) + interp = assemble(interpolate(function, V)) return _bezier_plot(interp, axes, complex_component=complex_component, **kwargs) y_vals = _bezier_calculate_points(function) diff --git a/firedrake/utility_meshes.py b/firedrake/utility_meshes.py index 6540611752..11bab9fc0f 100644 --- a/firedrake/utility_meshes.py +++ b/firedrake/utility_meshes.py @@ -11,7 +11,7 @@ Function, Constant, assemble, - Interpolate, + interpolate, FiniteElement, interval, tetrahedron, @@ -2351,7 +2351,7 @@ def OctahedralSphereMesh( ) if degree > 1: # use it to build a higher-order mesh - m = assemble(Interpolate(ufl.SpatialCoordinate(m), VectorFunctionSpace(m, "CG", degree))) + m = assemble(interpolate(ufl.SpatialCoordinate(m), VectorFunctionSpace(m, "CG", degree))) m = mesh.Mesh( m, name=name, @@ -2386,11 +2386,11 @@ def OctahedralSphereMesh( # Make a copy of the coordinates so that we can blend two different # mappings near the pole Vc = m.coordinates.function_space() - Xlatitudinal = assemble(Interpolate( + Xlatitudinal = assemble(interpolate( Constant(radius) * ufl.as_vector([x * scale, y * scale, znew]), Vc )) Vlow = VectorFunctionSpace(m, "CG", 1) - Xlow = assemble(Interpolate(Xlatitudinal, Vlow)) + Xlow = assemble(interpolate(Xlatitudinal, Vlow)) r = ufl.sqrt(Xlow[0] ** 2 + Xlow[1] ** 2 + Xlow[2] ** 2) Xradial = Constant(radius) * Xlow / r diff --git a/tests/firedrake/adjoint/test_reduced_functional.py b/tests/firedrake/adjoint/test_reduced_functional.py index acd07d259f..d4f9a863d0 100644 --- a/tests/firedrake/adjoint/test_reduced_functional.py +++ b/tests/firedrake/adjoint/test_reduced_functional.py @@ -214,7 +214,7 @@ def test_interpolate(): f = Function(V) f.dat.data[:] = 2 - J = assemble(Interpolate(f**2, c)) + J = assemble(interpolate(f**2, c)) Jhat = ReducedFunctional(J, Control(f)) h = Function(V) @@ -244,7 +244,7 @@ def test_interpolate_mixed(): f1, f2 = split(f) exprs = [f2 * div(f1)**2, grad(f2) * div(f1)] expr = as_vector([e[i] for e in exprs for i in np.ndindex(e.ufl_shape)]) - J = assemble(Interpolate(expr, c)) + J = assemble(interpolate(expr, c)) Jhat = ReducedFunctional(J, Control(f)) h = Function(V) diff --git a/tests/firedrake/external_operators/test_external_operators.py b/tests/firedrake/external_operators/test_external_operators.py index b6153f3d1f..a47953a566 100644 --- a/tests/firedrake/external_operators/test_external_operators.py +++ b/tests/firedrake/external_operators/test_external_operators.py @@ -104,7 +104,7 @@ def test_assemble(V, f): assert isinstance(jac, MatrixBase) # Assemble the exact Jacobian, i.e. the interpolation matrix: `Interpolate(dexpr(u,v,w)/du, V)` - jac_exact = assemble(Interpolate(derivative(expr(u, v, w), u), V)) + jac_exact = assemble(interpolate(derivative(expr(u, v, w), u), V)) np.allclose(jac.petscmat[:, :], jac_exact.petscmat[:, :], rtol=1e-14) # -- dNdu(u, v, w; δu, v*) (TLM) -- # diff --git a/tests/firedrake/multigrid/test_poisson_gtmg.py b/tests/firedrake/multigrid/test_poisson_gtmg.py index f70e5c6825..a4154dd392 100644 --- a/tests/firedrake/multigrid/test_poisson_gtmg.py +++ b/tests/firedrake/multigrid/test_poisson_gtmg.py @@ -60,7 +60,7 @@ def p1_callback(): if custom_transfer: P1 = get_p1_space() V = FunctionSpace(mesh, "DGT", degree - 1) - I = assemble(Interpolate(TrialFunction(P1), V)).petscmat + I = assemble(interpolate(TrialFunction(P1), V)).petscmat R = PETSc.Mat().createTranspose(I) appctx['interpolation_matrix'] = I appctx['restriction_matrix'] = R diff --git a/tests/firedrake/regression/test_adjoint_operators.py b/tests/firedrake/regression/test_adjoint_operators.py index 03557bf435..d0212636c4 100644 --- a/tests/firedrake/regression/test_adjoint_operators.py +++ b/tests/firedrake/regression/test_adjoint_operators.py @@ -729,7 +729,7 @@ def test_copy_function(): g = f.copy(deepcopy=True) J = assemble(g*dx) rf = ReducedFunctional(J, Control(f)) - a = assemble(Interpolate(-one, V)) + a = assemble(interpolate(-one, V)) assert np.isclose(rf(a), -J) diff --git a/tests/firedrake/regression/test_interpolate_zany.py b/tests/firedrake/regression/test_interpolate_zany.py index b2054843cd..875d0fb0c1 100644 --- a/tests/firedrake/regression/test_interpolate_zany.py +++ b/tests/firedrake/regression/test_interpolate_zany.py @@ -117,7 +117,7 @@ def test_interpolate_zany_into_vom(V, mesh, which, expr_at_vom): P0 = expr_at_vom.function_space() # Interpolate a Function into P0(vom) - f_at_vom = assemble(Interpolate(fexpr, P0)) + f_at_vom = assemble(interpolate(fexpr, P0)) assert numpy.allclose(f_at_vom.dat.data_ro, expr_at_vom.dat.data_ro) # Construct a Cofunction on P0(vom)* @@ -125,10 +125,10 @@ def test_interpolate_zany_into_vom(V, mesh, which, expr_at_vom): expected_action = assemble(action(Fvom, expr_at_vom)) # Interpolate a Function into Fvom - f_at_vom = assemble(Interpolate(fexpr, Fvom)) + f_at_vom = assemble(interpolate(fexpr, Fvom)) assert numpy.allclose(f_at_vom, expected_action) # Interpolate a TestFunction into Fvom - expr_vom = assemble(Interpolate(vexpr, Fvom)) + expr_vom = assemble(interpolate(vexpr, Fvom)) f_at_vom = assemble(action(expr_vom, f)) assert numpy.allclose(f_at_vom, expected_action) diff --git a/tests/firedrake/submesh/test_submesh_interpolate.py b/tests/firedrake/submesh/test_submesh_interpolate.py index a26c1acb08..92e422cf6a 100644 --- a/tests/firedrake/submesh/test_submesh_interpolate.py +++ b/tests/firedrake/submesh/test_submesh_interpolate.py @@ -51,7 +51,7 @@ def _test_submesh_interpolate_cell_cell(mesh, subdomain_cond, fe_fesub): f = Function(V_).interpolate(f) v0 = Coargument(V.dual(), 0) v1 = TrialFunction(Vsub) - interp = Interpolate(v1, v0, allow_missing_dofs=True) + interp = interpolate(v1, v0, allow_missing_dofs=True) A = assemble(interp) g = assemble(action(A, gsub)) assert assemble(inner(g - f, g - f) * dx(label_value)).real < 1e-14 @@ -165,7 +165,7 @@ def test_submesh_interpolate_subcell_subcell_2_processes(): f_l.dat.data_with_halos[:] = 3.0 v0 = Coargument(V_r.dual(), 0) v1 = TrialFunction(V_l) - interp = Interpolate(v1, v0, allow_missing_dofs=True) + interp = interpolate(v1, v0, allow_missing_dofs=True) A = assemble(interp) f_r = assemble(action(A, f_l)) g_r = Function(V_r).interpolate(conditional(x < 2.001, 3.0, 0.0)) From 0c4ce97aab6bb10940533e0a957671c52fe77421 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 23 Sep 2025 17:28:47 +0100 Subject: [PATCH 05/21] update `test_interp_dual.py` --- .../firedrake/regression/test_interp_dual.py | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index 50e29b05cb..ccd4de13f0 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -54,7 +54,7 @@ def test_assemble_interp_adjoint_tensor(mesh, V1, f1): def test_assemble_interp_operator(V2, f1): # Check type - If1 = Interpolate(f1, V2) + If1 = Interpolate(f1, Argument(V2.dual(), 0)) assert isinstance(If1, ufl.Interpolate) # -- I(f1, V2) -- # @@ -89,7 +89,7 @@ def test_assemble_interp_matrix(V1, V2, f1): def test_assemble_interp_tlm(V1, V2, f1): # -- Action(I(v1, V2), f1) -- # v1 = TrialFunction(V1) - Iv1 = Interpolate(v1, V2) + Iv1 = Interpolate(v1, Argument(V2.dual(), 0)) b = assemble(interpolate(f1, V2)) assembled_action_Iv1 = assemble(action(Iv1, f1)) @@ -99,7 +99,7 @@ def test_assemble_interp_tlm(V1, V2, f1): def test_assemble_interp_adjoint_matrix(V1, V2): # -- Adjoint(I(v1, V2)) -- # v1 = TrialFunction(V1) - Iv1 = Interpolate(v1, V2) + Iv1 = Interpolate(v1, Argument(V2.dual(), 0)) v2 = TestFunction(V2) c2 = assemble(conj(v2) * dx) @@ -120,7 +120,7 @@ def test_assemble_interp_adjoint_matrix(V1, V2): def test_assemble_interp_adjoint_model(V1, V2): # -- Action(Adjoint(I(v1, v2)), fstar) -- # v1 = TrialFunction(V1) - Iv1 = Interpolate(v1, V2) + Iv1 = Interpolate(v1, Argument(V2.dual(), 0)) fstar = Cofunction(V2.dual()) v = Argument(V1, 0) @@ -167,9 +167,9 @@ def test_assemble_base_form_operator_expressions(mesh): f2 = Function(V1).interpolate(sin(2*pi*y)) f3 = Function(V1).interpolate(cos(2*pi*x)) - If1 = Interpolate(f1, V2) - If2 = Interpolate(f2, V2) - If3 = Interpolate(f3, V2) + If1 = Interpolate(f1, Argument(V2.dual(), 0)) + If2 = Interpolate(f2, Argument(V2.dual(), 0)) + If3 = Interpolate(f3, Argument(V2.dual(), 0)) # Sum of BaseFormOperators (1-form) res = assemble(If1 + If2 + If3) @@ -234,7 +234,7 @@ def test_solve_interp_f(mesh): # -- Solution where the source term is interpolated via `ufl.Interpolate` u2 = Function(V1) - If = Interpolate(f1, V2) + If = Interpolate(f1, Argument(V2.dual(), 0)) # This requires assembling If F2 = inner(grad(u2), grad(w))*dx + inner(u2, w)*dx - inner(If, w)*dx solve(F2 == 0, u2) @@ -267,7 +267,7 @@ def test_solve_interp_u(mesh): # -- Solution where u2 is interpolated via `ufl.Interpolate` (mat-free) u2 = Function(V1) # Iu is the identity - Iu = Interpolate(u2, V1) + Iu = Interpolate(u2, Argument(V1.dual(), 0)) # This requires assembling the action the Jacobian of Iu F2 = inner(grad(u2), grad(w))*dx + inner(Iu, w)*dx - inner(f, w)*dx solve(F2 == 0, u2, solver_parameters={"mat_type": "matfree", @@ -278,7 +278,7 @@ def test_solve_interp_u(mesh): # Same problem with grad(Iu) instead of grad(Iu) u2 = Function(V1) # Iu is the identity - Iu = Interpolate(u2, V1) + Iu = Interpolate(u2, Argument(V1.dual(), 0)) # This requires assembling the action the Jacobian of Iu F2 = inner(grad(Iu), grad(w))*dx + inner(Iu, w)*dx - inner(f, w)*dx solve(F2 == 0, u2, solver_parameters={"mat_type": "matfree", From cf815919aabf12ca3177f2506743da1277c5bb09 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 23 Sep 2025 17:35:47 +0100 Subject: [PATCH 06/21] DROP BEFORE MERGE: use UFL branch --- .github/workflows/core.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index 55df1c8e01..61ea7f9d34 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -210,6 +210,7 @@ jobs: --extra-index-url https://download.pytorch.org/whl/cpu \ "$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]" + pip install -I "fenics-ufl @ git+https://github.com/firedrakeproject/ufl.git@leo/sanitise-interpolate" firedrake-clean pip list From c1b93f62332fb65f49e908815c97109318c69b9b Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 23 Sep 2025 17:36:30 +0100 Subject: [PATCH 07/21] move FunctionSpace check into `Interpolate` fix --- firedrake/interpolation.py | 27 +++++++-------------------- 1 file changed, 7 insertions(+), 20 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index bf3fcb5013..7c9258000b 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -48,7 +48,7 @@ class Interpolate(ufl.Interpolate): - def __init__(self, expr, v, + def __init__(self, expr, V, subset=None, access=None, allow_missing_dofs=False, @@ -60,7 +60,7 @@ def __init__(self, expr, v, ---------- expr : ufl.core.expr.Expr or ufl.BaseForm The UFL expression to interpolate. - v : firedrake.functionspaceimpl.WithGeometryBase or firedrake.ufl_expr.Coargument + V : firedrake.functionspaceimpl.WithGeometryBase or firedrake.ufl_expr.Coargument The function space to interpolate into or the coargument defined on the dual of the function space to interpolate into. subset : pyop2.types.set.Subset @@ -95,20 +95,12 @@ def __init__(self, expr, v, between a VOM and its input ordering. Defaults to ``True`` which uses SF broadcast and reduce operations. """ - # Check function space - expr = ufl.as_ufl(expr) - if isinstance(v, functionspaceimpl.WithGeometry): - expr_args = extract_arguments(expr) + if isinstance(V, functionspaceimpl.WithGeometry): + # Need to create a Firedrake Argument so that it has a .function_space() method + expr_args = extract_arguments(ufl.as_ufl(expr)) is_adjoint = len(expr_args) and expr_args[0].number() == 0 - v = Argument(v.dual(), 1 if is_adjoint else 0) - - V = v.arguments()[0].function_space() - if len(expr.ufl_shape) != len(V.value_shape): - raise RuntimeError(f'Rank mismatch: Expression rank {len(expr.ufl_shape)}, FunctionSpace rank {len(V.value_shape)}') - - if expr.ufl_shape != V.value_shape: - raise RuntimeError('Shape mismatch: Expression shape {expr.ufl_shape}, FunctionSpace shape {V.value_shape}') - super().__init__(expr, v) + V = Argument(V.dual(), 1 if is_adjoint else 0) + super().__init__(expr, V) # -- Interpolate data (e.g. `subset` or `access`) -- # self.interp_data = {"subset": subset, @@ -174,11 +166,6 @@ def interpolate(expr, V, subset=None, access=None, allow_missing_dofs=False, def reduction (hence using MIN will compute the MIN between the existing values and any new values). """ - if isinstance(V, functionspaceimpl.WithGeometry): - # Need to create a Firedrake Argument so that it has a .function_space() method - expr_args = extract_arguments(ufl.as_ufl(expr)) - is_adjoint = len(expr_args) and expr_args[0].number() == 0 - V = Argument(V.dual(), 1 if is_adjoint else 0) return Interpolate( expr, V, subset=subset, access=access, allow_missing_dofs=allow_missing_dofs, default_missing_val=default_missing_val, matfree=matfree From 8f764c786d32491d8108a9ce4c97f3282a4826e2 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 23 Sep 2025 17:45:44 +0100 Subject: [PATCH 08/21] lint --- firedrake/interpolation.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 7c9258000b..70642cb105 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -27,11 +27,10 @@ import firedrake import firedrake.bcs from firedrake import tsfc_interface, utils, functionspaceimpl -from firedrake.ufl_expr import Argument, Coargument, action, adjoint as expr_adjoint +from firedrake.ufl_expr import Argument, action, adjoint as expr_adjoint 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 From e9e92dd66f0c9d21a12d35540691c87c729aefd6 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 24 Sep 2025 11:52:24 +0100 Subject: [PATCH 09/21] test -> trial test -> trial fix hypre-ads --- firedrake/preconditioners/hypre_ads.py | 6 +++--- firedrake/preconditioners/hypre_ams.py | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/firedrake/preconditioners/hypre_ads.py b/firedrake/preconditioners/hypre_ads.py index 89c10dc438..98443f2c75 100644 --- a/firedrake/preconditioners/hypre_ads.py +++ b/firedrake/preconditioners/hypre_ads.py @@ -1,7 +1,7 @@ from firedrake.preconditioners.base import PCBase from firedrake.petsc import PETSc from firedrake.function import Function -from firedrake.ufl_expr import TestFunction +from firedrake.ufl_expr import TrialFunction from firedrake.dmhooks import get_function_space from firedrake.preconditioners.hypre_ams import chop from firedrake.interpolation import interpolate @@ -31,12 +31,12 @@ def initialize(self, obj): NC1 = V.reconstruct(family="N1curl" if mesh.ufl_cell().is_simplex() else "NCE", degree=1) G_callback = appctx.get("get_gradient", None) if G_callback is None: - G = chop(assemble(interpolate(grad(TestFunction(P1)), NC1)).petscmat) + G = chop(assemble(interpolate(grad(TrialFunction(P1)), NC1)).petscmat) else: G = G_callback(P1, NC1) C_callback = appctx.get("get_curl", None) if C_callback is None: - C = chop(assemble(interpolate(curl(TestFunction(NC1)), V)).petscmat) + C = chop(assemble(interpolate(curl(TrialFunction(NC1)), V)).petscmat) else: C = C_callback(NC1, V) diff --git a/firedrake/preconditioners/hypre_ams.py b/firedrake/preconditioners/hypre_ams.py index 9a59702af4..594fe88590 100644 --- a/firedrake/preconditioners/hypre_ams.py +++ b/firedrake/preconditioners/hypre_ams.py @@ -2,7 +2,7 @@ from firedrake.preconditioners.base import PCBase from firedrake.petsc import PETSc from firedrake.function import Function -from firedrake.ufl_expr import TestFunction +from firedrake.ufl_expr import TrialFunction from firedrake.dmhooks import get_function_space from firedrake.utils import complex_mode from firedrake.interpolation import interpolate @@ -51,7 +51,7 @@ def initialize(self, obj): P1 = V.reconstruct(family="Lagrange", degree=1) G_callback = appctx.get("get_gradient", None) if G_callback is None: - G = chop(assemble(interpolate(grad(TestFunction(P1)), V)).petscmat) + G = chop(assemble(interpolate(grad(TrialFunction(P1)), V)).petscmat) else: G = G_callback(P1, V) From ba75cd5f0bba911b325d23c0d41823aa17beba4f Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 1 Oct 2025 13:53:33 +0100 Subject: [PATCH 10/21] tidy function.interpolate --- firedrake/function.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/firedrake/function.py b/firedrake/function.py index fed408c8cc..68290b8ce4 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -382,12 +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 = interpolate( - expression, V, subset=subset, allow_missing_dofs=allow_missing_dofs, - default_missing_val=default_missing_val - ) + interp = interpolate(expression, V, **kwargs) return assemble(interp, tensor=self, ad_block_tag=ad_block_tag) def zero(self, subset=None): From 83ef532ea6f2b58e1df803450746bf61f9c1fc72 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 1 Oct 2025 13:54:27 +0100 Subject: [PATCH 11/21] tidy cofunction.interpolate cofunction docstring --- firedrake/cofunction.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/firedrake/cofunction.py b/firedrake/cofunction.py index f0fda10f63..9ee867c622 100644 --- a/firedrake/cofunction.py +++ b/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 From 8adf841fb34f6b7112ae46c4cc790683ae7d2387 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 1 Oct 2025 14:10:49 +0100 Subject: [PATCH 12/21] tidy type hints in function.py --- firedrake/function.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/firedrake/function.py b/firedrake/function.py index 68290b8ce4..04579850a5 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -14,7 +14,7 @@ from numbers import Number from pathlib import Path from functools import partial -from typing import Tuple +from typing import Tuple, Optional from pyop2 import op2, mpi from pyop2.exceptions import DataTypeError, DataValueError @@ -362,7 +362,7 @@ def function_space(self): @PETSc.Log.EventDecorator() def interpolate(self, expression: ufl.classes.Expr, - ad_block_tag: str | None = None, + ad_block_tag: Optional[str] = None, **kwargs): """Interpolate an expression onto this :class:`Function`. @@ -697,13 +697,13 @@ 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: r"""Convenience class for evaluating a :class:`Function` at a set of points.""" - def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: float | None = None, + def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: Optional[float] = None, missing_points_behaviour: str = "error", redundant: bool = True) -> None: r""" Parameters @@ -712,7 +712,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 From 9f01a4b3cd5a788d20e3c6da5dcde7b46c91ea41 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 1 Oct 2025 15:04:55 +0100 Subject: [PATCH 13/21] remove UFL branch --- .github/workflows/core.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index 61ea7f9d34..55df1c8e01 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -210,7 +210,6 @@ jobs: --extra-index-url https://download.pytorch.org/whl/cpu \ "$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]" - pip install -I "fenics-ufl @ git+https://github.com/firedrakeproject/ufl.git@leo/sanitise-interpolate" firedrake-clean pip list From d8093ecd3e31582c45692c604762be8f819477b8 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 1 Oct 2025 15:19:01 +0100 Subject: [PATCH 14/21] lint --- firedrake/function.py | 6 +++--- firedrake/interpolation.py | 9 ++++----- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/firedrake/function.py b/firedrake/function.py index 04579850a5..38d6b6e76b 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -14,7 +14,7 @@ from numbers import Number from pathlib import Path from functools import partial -from typing import Tuple, Optional +from typing import Tuple from pyop2 import op2, mpi from pyop2.exceptions import DataTypeError, DataValueError @@ -362,7 +362,7 @@ def function_space(self): @PETSc.Log.EventDecorator() def interpolate(self, expression: ufl.classes.Expr, - ad_block_tag: Optional[str] = None, + ad_block_tag: str | None = None, **kwargs): """Interpolate an expression onto this :class:`Function`. @@ -703,7 +703,7 @@ def __str__(self): class PointEvaluator: r"""Convenience class for evaluating a :class:`Function` at a set of points.""" - def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: Optional[float] = None, + def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: float | None = None, missing_points_behaviour: str = "error", redundant: bool = True) -> None: r""" Parameters diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 70642cb105..c4c2ed238c 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -4,14 +4,13 @@ import abc import warnings from collections.abc import Iterable -from typing import Literal from functools import partial, singledispatch from typing import Hashable import FIAT import ufl import finat.ufl -from ufl.algorithms import extract_arguments, extract_coefficients, replace +from ufl.algorithms import extract_arguments, extract_coefficients from ufl.domain import as_domain, extract_unique_domain from pyop2 import op2 @@ -25,9 +24,9 @@ import finat import firedrake -import firedrake.bcs from firedrake import tsfc_interface, utils, functionspaceimpl -from firedrake.ufl_expr import Argument, action, adjoint as expr_adjoint +from firedrake.ufl_expr import Argument, Coargument, action, adjoint as expr_adjoint +from firedrake.cofunction import Cofunction from firedrake.mesh import MissingPointsBehaviour, VertexOnlyMeshMissingPointsError, VertexOnlyMeshTopology from firedrake.petsc import PETSc from firedrake.halo import _get_mtype as get_dat_mpi_type @@ -263,7 +262,7 @@ def __init__( V: ufl.FunctionSpace | firedrake.function.Function, subset: op2.Subset | None = None, freeze_expr: bool = False, - access: Literal[op2.WRITE, op2.MIN, op2.MAX, op2.INC] | None = None, + access: op2.Access | None = None, bcs: Iterable[firedrake.bcs.BCBase] | None = None, allow_missing_dofs: bool = False, matfree: bool = True From d15da7556c48aee16810f3c88246478e1069d3b2 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 1 Oct 2025 15:23:06 +0100 Subject: [PATCH 15/21] fix typing fix --- firedrake/interpolation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index c4c2ed238c..989f081dac 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -5,7 +5,7 @@ import warnings from collections.abc import Iterable from functools import partial, singledispatch -from typing import Hashable +from typing import Hashable, Literal import FIAT import ufl @@ -262,7 +262,7 @@ def __init__( V: ufl.FunctionSpace | firedrake.function.Function, subset: op2.Subset | None = None, freeze_expr: bool = False, - access: op2.Access | None = None, + access: Literal[op2.WRITE, op2.MIN, op2.MAX, op2.INC] | None = None, bcs: Iterable[firedrake.bcs.BCBase] | None = None, allow_missing_dofs: bool = False, matfree: bool = True From e9d6ba93f31b0164d3a14b4ec7e353bf294787ad Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 1 Oct 2025 22:43:50 +0100 Subject: [PATCH 16/21] runtimeerror -> valueerror --- tests/firedrake/regression/test_function.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/firedrake/regression/test_function.py b/tests/firedrake/regression/test_function.py index 85c6cdc0a0..1f9ed7a429 100644 --- a/tests/firedrake/regression/test_function.py +++ b/tests/firedrake/regression/test_function.py @@ -81,22 +81,22 @@ def test_firedrake_tensor_function_nonstandard_shape(W_nonstandard_shape): def test_mismatching_rank_interpolation(V): f = Function(V) - with pytest.raises(RuntimeError): + with pytest.raises(ValueError): f.interpolate(Constant((1, 2))) VV = VectorFunctionSpace(V.mesh(), 'CG', 1) f = Function(VV) - with pytest.raises(RuntimeError): + with pytest.raises(ValueError): f.interpolate(Constant((1, 2))) VVV = TensorFunctionSpace(V.mesh(), 'CG', 1) f = Function(VVV) - with pytest.raises(RuntimeError): + with pytest.raises(ValueError): f.interpolate(Constant((1, 2))) def test_mismatching_shape_interpolation(V): VV = VectorFunctionSpace(V.mesh(), 'CG', 1) f = Function(VV) - with pytest.raises(RuntimeError): + with pytest.raises(ValueError): f.interpolate(Constant([1] * (VV.value_shape[0] + 1))) From 52c6203c6c61e3fba9eccf7d2a8a6898f20e9624 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 1 Oct 2025 22:48:12 +0100 Subject: [PATCH 17/21] add check for shape mismatch to `Interpolate` --- firedrake/interpolation.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 989f081dac..9341f1dd8d 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -98,6 +98,11 @@ def __init__(self, expr, V, expr_args = extract_arguments(ufl.as_ufl(expr)) is_adjoint = len(expr_args) and expr_args[0].number() == 0 V = Argument(V.dual(), 1 if is_adjoint else 0) + + target_shape = V.arguments()[0].function_space().value_shape + if expr.ufl_shape != target_shape: + raise ValueError(f"Shape mismatch: Expression shape {expr.ufl_shape}, FunctionSpace shape {target_shape}.") + super().__init__(expr, V) # -- Interpolate data (e.g. `subset` or `access`) -- # From aff95243cb1453a5b0793252a6b627bb56d7fe66 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Thu, 2 Oct 2025 11:24:03 +0100 Subject: [PATCH 18/21] use ufl.as_expr --- firedrake/interpolation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 9341f1dd8d..dee24dc64c 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -93,9 +93,10 @@ def __init__(self, expr, V, between a VOM and its input ordering. Defaults to ``True`` which uses SF broadcast and reduce operations. """ + expr = ufl.as_ufl(expr) if isinstance(V, functionspaceimpl.WithGeometry): # Need to create a Firedrake Argument so that it has a .function_space() method - expr_args = extract_arguments(ufl.as_ufl(expr)) + expr_args = extract_arguments(expr) is_adjoint = len(expr_args) and expr_args[0].number() == 0 V = Argument(V.dual(), 1 if is_adjoint else 0) From 93406b85486abc3f7b0a943fc25b0ecb0a9ebe43 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 6 Oct 2025 12:34:27 +0100 Subject: [PATCH 19/21] update expr arg check --- firedrake/interpolation.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index dee24dc64c..ab908d2a01 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -12,6 +12,7 @@ import finat.ufl from ufl.algorithms import extract_arguments, extract_coefficients from ufl.domain import as_domain, extract_unique_domain +from ufl.duals import is_dual from pyop2 import op2 from pyop2.caching import memory_and_disk_cache @@ -96,8 +97,8 @@ def __init__(self, expr, V, expr = ufl.as_ufl(expr) if isinstance(V, functionspaceimpl.WithGeometry): # Need to create a Firedrake Argument so that it has a .function_space() method - expr_args = extract_arguments(expr) - is_adjoint = len(expr_args) and expr_args[0].number() == 0 + expr_arg_numbers = {arg.number() for arg in extract_arguments(expr) if not is_dual(arg)} + is_adjoint = len(expr_arg_numbers) and expr_arg_numbers[0] == 0 V = Argument(V.dual(), 1 if is_adjoint else 0) target_shape = V.arguments()[0].function_space().value_shape From efedc48aaf3c5ce2616bc0329c2f8bf2dbb4fb06 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 6 Oct 2025 12:35:15 +0100 Subject: [PATCH 20/21] lint --- firedrake/interpolation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index ab908d2a01..d58323393a 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -27,7 +27,6 @@ import firedrake from firedrake import tsfc_interface, utils, functionspaceimpl from firedrake.ufl_expr import Argument, Coargument, action, adjoint as expr_adjoint -from firedrake.cofunction import Cofunction from firedrake.mesh import MissingPointsBehaviour, VertexOnlyMeshMissingPointsError, VertexOnlyMeshTopology from firedrake.petsc import PETSc from firedrake.halo import _get_mtype as get_dat_mpi_type From f8c2318b9ae19bdb6d1a28c7e2729ca675c62571 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 6 Oct 2025 12:45:19 +0100 Subject: [PATCH 21/21] fix check --- firedrake/interpolation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index d58323393a..fb60ed7da2 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -97,7 +97,7 @@ def __init__(self, expr, V, if isinstance(V, functionspaceimpl.WithGeometry): # Need to create a Firedrake Argument so that it has a .function_space() method expr_arg_numbers = {arg.number() for arg in extract_arguments(expr) if not is_dual(arg)} - is_adjoint = len(expr_arg_numbers) and expr_arg_numbers[0] == 0 + is_adjoint = len(expr_arg_numbers) and expr_arg_numbers == {0} V = Argument(V.dual(), 1 if is_adjoint else 0) target_shape = V.arguments()[0].function_space().value_shape