From 57784e152ce3a96b60923239706c1d0e096fc7a7 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Fri, 19 Dec 2025 16:35:20 +0000 Subject: [PATCH 01/34] enable tests --- .../regression/test_interpolate_cross_mesh.py | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index 2011dc9466..0558a2f318 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -49,14 +49,7 @@ def make_high_order(m_low_order, degree): "unitsquare_vfs", "unitsquare_tfs", "unitsquare_N1curl_source", - pytest.param( - "unitsquare_SminusDiv_destination", - marks=pytest.mark.xfail( - # CalledProcessError is so the parallel tests correctly xfail - raises=(subprocess.CalledProcessError, NotImplementedError), - reason="Can only interpolate into spaces with point evaluation nodes", - ), - ), + "unitsquare_SminusDiv_destination", "unitsquare_Regge_source", # This test fails in complex mode pytest.param("spheresphere", marks=pytest.mark.skipcomplex), @@ -438,10 +431,9 @@ def test_interpolate_cross_mesh_not_point_eval(): atol = 1e-8 # default # This might not make much mathematical sense, but it should test if we get # the not implemented error for non-point evaluation nodes! - with pytest.raises(NotImplementedError): - interpolate_function( - m_src, m_dest, V_src, V_dest, dest_eval, expected, expr_src, expr_dest, atol - ) + interpolate_function( + m_src, m_dest, V_src, V_dest, dest_eval, expected, expr_src, expr_dest, atol + ) def interpolate_function( From 38ece6b3113b5274766cf42fa3089bb8b8b8d4b1 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Sun, 21 Dec 2025 18:30:05 +0000 Subject: [PATCH 02/34] remove check --- firedrake/interpolation.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 6f0954d8ea..b69fc22647 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -422,17 +422,6 @@ def __init__(self, expr: Interpolate): else: self.access = op2.WRITE - # TODO check V.finat_element.is_lagrange() once https://github.com/firedrakeproject/fiat/pull/200 is released - target_element = self.target_space.ufl_element() - if not ((isinstance(target_element, MixedElement) - and all(sub.mapping() == "identity" for sub in target_element.sub_elements)) - or target_element.mapping() == "identity"): - # Identity mapping between reference cell and physical coordinates - # implies point evaluation nodes. - raise NotImplementedError( - "Can only cross-mesh interpolate into spaces with point evaluation nodes." - ) - if self.allow_missing_dofs: self.missing_points_behaviour = MissingPointsBehaviour.IGNORE else: From 013f31a2b6a3192b768b02f08f98a7c16e65c868 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Sun, 21 Dec 2025 19:27:35 +0000 Subject: [PATCH 03/34] add test --- .../test_cross_mesh_non_lagrange.py | 78 +++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 tests/firedrake/regression/test_cross_mesh_non_lagrange.py diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py new file mode 100644 index 0000000000..2540930274 --- /dev/null +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -0,0 +1,78 @@ +from firedrake import * +import pytest +import numpy as np +from finat.quadrature import QuadratureRule + +def make_quadrature_space(V): + """Builds a quadrature space on the target mesh. + + This space has point evaluation dofs at the quadrature PointSet + of the target space's element + + """ + fe = V.finat_element + _, ps = fe.dual_basis + wts = np.full(len(ps.points), np.nan) # These can be any number since we never integrate + scheme = QuadratureRule(ps, wts, ref_el=fe.cell) + element = FiniteElement("Quadrature", degree=fe.degree, quad_scheme=scheme) + return VectorFunctionSpace(V.mesh(), element) + + +@pytest.fixture +def mesh1(): + return UnitSquareMesh(5, 5) + +@pytest.fixture +def mesh2(): + return UnitSquareMesh(3, 3) + +@pytest.mark.parametrize("element,degree", [("RT", 1), ("RT", 2), ("RT", 3), + ("BDM", 1), ("BDM", 2), ("BDM", 3), + ("BDFM", 2)]) +def test_hdiv_cross_mesh_oneform(mesh1, mesh2, element, degree): + V_source = VectorFunctionSpace(mesh1, "CG", 2) + V_target = FunctionSpace(mesh2, element, degree) + + x, y = SpatialCoordinate(mesh1) + f_source = Function(V_source).interpolate(as_vector([x, y])) + + # Make intermediate Quadrature space on target mesh + Q_target = make_quadrature_space(V_target) + + # Interp V_source -> Q + I1 = interpolate(f_source, Q_target) + f_quadrature = assemble(I1) + + # Interp Q -> V_target + I2 = interpolate(f_quadrature, V_target) + f_target = assemble(I2) + + x1, y1 = SpatialCoordinate(mesh2) + f_direct = Function(V_target).interpolate(as_vector([x1, y1])) + + assert np.allclose(f_target.dat.data_ro, f_direct.dat.data_ro) + +@pytest.mark.parametrize("element,degree", [("N1curl", 1), ("N1curl", 2), ("N1curl", 3), + ("N2curl", 1), ("N2curl", 2), ("N2curl", 3),]) +def test_hcurl_cross_mesh_oneform(mesh1, mesh2, element, degree): + V_source = VectorFunctionSpace(mesh1, "CG", 2) + V_target = FunctionSpace(mesh2, element, degree) + + x, y = SpatialCoordinate(mesh1) + f_source = Function(V_source).interpolate(as_vector([x, y])) + + # Make intermediate Quadrature space on target mesh + Q_target = make_quadrature_space(V_target) + + # Interp V_source -> Q + I1 = interpolate(f_source, Q_target) + f_quadrature = assemble(I1) + + # Interp Q -> V_target + I2 = interpolate(f_quadrature, V_target) + f_target = assemble(I2) + + x1, y1 = SpatialCoordinate(mesh2) + f_direct = Function(V_target).interpolate(as_vector([x1, y1])) + + assert np.allclose(f_target.dat.data_ro, f_direct.dat.data_ro) From ee75c503798de3fc893670e789dab825ce416d0c Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 22 Dec 2025 17:13:55 +0000 Subject: [PATCH 04/34] update test --- .../test_cross_mesh_non_lagrange.py | 80 +++++++++---------- 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 2540930274..6c88cf57f9 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -2,6 +2,18 @@ import pytest import numpy as np from finat.quadrature import QuadratureRule +from functools import partial + + +def fs_shape(V): + shape = V.ufl_function_space().value_shape + if len(shape) == 1: + fs_type = partial(VectorFunctionSpace, dim=shape[0]) + elif len(shape) == 2: + fs_type = partial(TensorFunctionSpace, shape=shape) + else: + raise ValueError("Invalid function space shape") + return fs_type def make_quadrature_space(V): """Builds a quadrature space on the target mesh. @@ -15,51 +27,40 @@ def make_quadrature_space(V): wts = np.full(len(ps.points), np.nan) # These can be any number since we never integrate scheme = QuadratureRule(ps, wts, ref_el=fe.cell) element = FiniteElement("Quadrature", degree=fe.degree, quad_scheme=scheme) - return VectorFunctionSpace(V.mesh(), element) - + return fs_shape(V)(V.mesh(), element) -@pytest.fixture -def mesh1(): - return UnitSquareMesh(5, 5) -@pytest.fixture -def mesh2(): - return UnitSquareMesh(3, 3) +@pytest.fixture(params=[("RT", 1), ("RT", 2), ("RT", 3), ("BDM", 1), ("BDM", 2), ("BDM", 3), + ("BDFM", 2), ("HHJ", 2),("N1curl", 1), ("N1curl", 2), ("N1curl", 3), + ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 1), ("GLS", 2), + ("GLS", 3), ("GLS2", 1), ("GLS2", 2), ("GLS2", 3)], + ids=lambda x: f"{x[0]}_{x[1]}") +def V_target(request): + element, degree = request.param + mesh = UnitSquareMesh(3, 3) + return FunctionSpace(mesh, element, degree) -@pytest.mark.parametrize("element,degree", [("RT", 1), ("RT", 2), ("RT", 3), - ("BDM", 1), ("BDM", 2), ("BDM", 3), - ("BDFM", 2)]) -def test_hdiv_cross_mesh_oneform(mesh1, mesh2, element, degree): - V_source = VectorFunctionSpace(mesh1, "CG", 2) - V_target = FunctionSpace(mesh2, element, degree) +def test_cross_mesh_oneform(V_target): + mesh1 = UnitSquareMesh(5, 5) + mesh2 = V_target.mesh() x, y = SpatialCoordinate(mesh1) - f_source = Function(V_source).interpolate(as_vector([x, y])) - - # Make intermediate Quadrature space on target mesh - Q_target = make_quadrature_space(V_target) - - # Interp V_source -> Q - I1 = interpolate(f_source, Q_target) - f_quadrature = assemble(I1) - - # Interp Q -> V_target - I2 = interpolate(f_quadrature, V_target) - f_target = assemble(I2) - x1, y1 = SpatialCoordinate(mesh2) - f_direct = Function(V_target).interpolate(as_vector([x1, y1])) - - assert np.allclose(f_target.dat.data_ro, f_direct.dat.data_ro) - -@pytest.mark.parametrize("element,degree", [("N1curl", 1), ("N1curl", 2), ("N1curl", 3), - ("N2curl", 1), ("N2curl", 2), ("N2curl", 3),]) -def test_hcurl_cross_mesh_oneform(mesh1, mesh2, element, degree): - V_source = VectorFunctionSpace(mesh1, "CG", 2) - V_target = FunctionSpace(mesh2, element, degree) - x, y = SpatialCoordinate(mesh1) - f_source = Function(V_source).interpolate(as_vector([x, y])) + shape = V_target.ufl_function_space().value_shape + if len(shape) == 1: + fs_type = partial(VectorFunctionSpace, dim=shape[0]) + expr1 = as_vector([x, y]) + expr2 = as_vector([x1, y1]) + elif len(shape) == 2: + fs_type = partial(TensorFunctionSpace, shape=shape) + expr1 = as_tensor([[x, x*y], [x*y, y]]) + expr2 = as_tensor([[x1, x1*y1], [x1*y1, y1]]) + else: + raise ValueError("Unsupported target space shape") + + V_source = fs_type(mesh1, "CG", 2) + f_source = Function(V_source).interpolate(expr1) # Make intermediate Quadrature space on target mesh Q_target = make_quadrature_space(V_target) @@ -72,7 +73,6 @@ def test_hcurl_cross_mesh_oneform(mesh1, mesh2, element, degree): I2 = interpolate(f_quadrature, V_target) f_target = assemble(I2) - x1, y1 = SpatialCoordinate(mesh2) - f_direct = Function(V_target).interpolate(as_vector([x1, y1])) + f_direct = Function(V_target).interpolate(expr2) assert np.allclose(f_target.dat.data_ro, f_direct.dat.data_ro) From c26ae0b19dfd4a0e482e85282efa12da0dee4fe4 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 5 Jan 2026 17:56:10 +0000 Subject: [PATCH 05/34] update tests --- .../test_cross_mesh_non_lagrange.py | 48 ++++++++++++++++--- .../regression/test_interpolate_cross_mesh.py | 4 +- 2 files changed, 43 insertions(+), 9 deletions(-) diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 6c88cf57f9..65032bdddb 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -35,19 +35,21 @@ def make_quadrature_space(V): ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 1), ("GLS", 2), ("GLS", 3), ("GLS2", 1), ("GLS2", 2), ("GLS2", 3)], ids=lambda x: f"{x[0]}_{x[1]}") -def V_target(request): +def V(request): element, degree = request.param mesh = UnitSquareMesh(3, 3) return FunctionSpace(mesh, element, degree) +# V_source -> Q -> V_target +# V_target^* -> Q^* -> V_source^* -def test_cross_mesh_oneform(V_target): +def test_cross_mesh_oneform(V): mesh1 = UnitSquareMesh(5, 5) - mesh2 = V_target.mesh() + mesh2 = V.mesh() x, y = SpatialCoordinate(mesh1) x1, y1 = SpatialCoordinate(mesh2) - shape = V_target.ufl_function_space().value_shape + shape = V.ufl_function_space().value_shape if len(shape) == 1: fs_type = partial(VectorFunctionSpace, dim=shape[0]) expr1 = as_vector([x, y]) @@ -63,16 +65,48 @@ def test_cross_mesh_oneform(V_target): f_source = Function(V_source).interpolate(expr1) # Make intermediate Quadrature space on target mesh - Q_target = make_quadrature_space(V_target) + Q_target = make_quadrature_space(V) # Interp V_source -> Q I1 = interpolate(f_source, Q_target) f_quadrature = assemble(I1) # Interp Q -> V_target - I2 = interpolate(f_quadrature, V_target) + I2 = interpolate(f_quadrature, V) f_target = assemble(I2) - f_direct = Function(V_target).interpolate(expr2) + f_direct = Function(V).interpolate(expr2) assert np.allclose(f_target.dat.data_ro, f_direct.dat.data_ro) + + +def test_cross_mesh_oneform_adjoint(V): + # Can already do Lagrange -> RT adjoint + # V^* -> Q^* -> V_target^* + mesh1 = UnitSquareMesh(7, 7) + x1 = SpatialCoordinate(mesh1) + V_target = fs_shape(V)(mesh1, "CG", 2) + + mesh2 = V.mesh() + x2 = SpatialCoordinate(mesh2) + + oneform_V = inner(x2, TestFunction(V)) * dx + + Q_target = make_quadrature_space(V) + + # Interp V^* -> Q^* + I1_adj = interpolate(TestFunction(Q_target), oneform_V) + cofunc_Q = assemble(I1_adj) + + # Interp Q^* -> V_target^* + I2_adj = interpolate(TestFunction(V_target), cofunc_Q) + cofunc_V = assemble(I2_adj) + + # cofunc_V = assemble(interpolate(TestFunction(V_target), oneform_target)) # V^* -> V_target^* + + cofunc_V_direct = assemble(inner(x1, TestFunction(V_target)) * dx) + + assert np.allclose(cofunc_V.dat.data_ro, cofunc_V_direct.dat.data_ro) + +if __name__ == "__main__": + pytest.main([__file__ + "::test_cross_mesh_oneform_adjoint[RT_1]"]) \ No newline at end of file diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index 0558a2f318..c57a8e025e 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -187,8 +187,8 @@ def parameters(request): expr_dest = 2 * SpatialCoordinate(m_dest) expected = 2 * coords V_src = VectorFunctionSpace(m_src, "CG", 2) - V_dest = FunctionSpace(m_dest, "SminusDiv", 2) # Not point evaluation nodes - V_dest_2 = FunctionSpace(m_dest, "SminusCurl", 2) # Not point evaluation nodes + V_dest = FunctionSpace(m_dest, "RT", 2) # Not point evaluation nodes + V_dest_2 = FunctionSpace(m_dest, "N1Curl", 2) # Not point evaluation nodes elif request.param == "unitsquare_Regge_source": m_src, m_dest, coords = unitsquaresetup() expr_src = outer(SpatialCoordinate(m_src), SpatialCoordinate(m_src)) From 4374ac451eba3f0c729ae6f34cf658f3c187afe5 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 5 Jan 2026 17:56:28 +0000 Subject: [PATCH 06/34] WIP need to map replaced coefficient correctly --- tsfc/driver.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/tsfc/driver.py b/tsfc/driver.py index f83db1d944..d6528b779f 100644 --- a/tsfc/driver.py +++ b/tsfc/driver.py @@ -245,15 +245,20 @@ def compile_expression_dual_evaluation(expression, ufl_element, *, complex_mode = is_complex(parameters["scalar_type"]) orig_coefficients = extract_coefficients(expression) - if isinstance(expression, ufl.Interpolate): - v, operand = expression.argument_slots() - else: - operand = expression - v = ufl.FunctionSpace(extract_unique_domain(operand), ufl_element) + v, operand = expression.argument_slots() # Map into reference space operand = apply_mapping(operand, ufl_element, domain) + if ufl_element.mapping() != "identity": + # Need to map dual argument for adjoint interpolation + ref_element = finat.ufl.WithMapping(ufl_element, "identity") + V = ufl.FunctionSpace(domain, ref_element) + if isinstance(v, ufl.Coargument): + v = ufl.Coargument(V.dual(), v.number()) + else: + v = ufl.Cofunction(V.dual()) + # Apply UFL preprocessing operand = ufl_utils.preprocess_expression(operand, complex_mode=complex_mode) operand = simplify_abs(operand, complex_mode) @@ -273,9 +278,6 @@ def compile_expression_dual_evaluation(expression, ufl_element, *, assert len(argument_multiindices) == len(arguments) # Replace coordinates (if any) unless otherwise specified by kwarg - if domain is None: - domain = extract_unique_domain(expression) - assert domain is not None builder._domain_integral_type_map = {domain: "cell"} builder._entity_ids = {domain: (0,)} From 37f023661129713f952d914adb5bcbcd2bc4dc7a Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 20 Jan 2026 10:10:50 +0000 Subject: [PATCH 07/34] fix test fixes --- .../test_cross_mesh_non_lagrange.py | 33 ++++++++++++------- tsfc/driver.py | 18 +++++----- 2 files changed, 30 insertions(+), 21 deletions(-) diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 65032bdddb..6651eada43 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -4,6 +4,8 @@ from finat.quadrature import QuadratureRule from functools import partial +from ufl.compound_expressions import deviatoric_expr_2x2 + def fs_shape(V): shape = V.ufl_function_space().value_shape @@ -30,14 +32,14 @@ def make_quadrature_space(V): return fs_shape(V)(V.mesh(), element) -@pytest.fixture(params=[("RT", 1), ("RT", 2), ("RT", 3), ("BDM", 1), ("BDM", 2), ("BDM", 3), - ("BDFM", 2), ("HHJ", 2),("N1curl", 1), ("N1curl", 2), ("N1curl", 3), - ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 1), ("GLS", 2), - ("GLS", 3), ("GLS2", 1), ("GLS2", 2), ("GLS2", 3)], +@pytest.fixture(params=[("RT", 2), ("RT", 3), ("RT", 4), ("BDM", 1), ("BDM", 2), ("BDM", 3), + ("BDFM", 2), ("HHJ", 2),("N1curl", 2), ("N1curl", 3), ("N1curl", 4), + ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 2), ("GLS", 3), + ("GLS", 4), ("GLS2", 1), ("GLS2", 2), ("GLS2", 3)], ids=lambda x: f"{x[0]}_{x[1]}") def V(request): element, degree = request.param - mesh = UnitSquareMesh(3, 3) + mesh = UnitSquareMesh(8, 8) return FunctionSpace(mesh, element, degree) # V_source -> Q -> V_target @@ -83,14 +85,24 @@ def test_cross_mesh_oneform(V): def test_cross_mesh_oneform_adjoint(V): # Can already do Lagrange -> RT adjoint # V^* -> Q^* -> V_target^* - mesh1 = UnitSquareMesh(7, 7) + mesh1 = UnitSquareMesh(2, 2) x1 = SpatialCoordinate(mesh1) - V_target = fs_shape(V)(mesh1, "CG", 2) + V_target = fs_shape(V)(mesh1, "CG", 1) mesh2 = V.mesh() x2 = SpatialCoordinate(mesh2) - oneform_V = inner(x2, TestFunction(V)) * dx + if len(V.value_shape) > 1: + expr = outer(x2, x2) + target_expr = outer(x1, x1) + if V.ufl_element().mapping() == "covariant contravariant Piola": + expr = deviatoric_expr_2x2(expr) + target_expr = deviatoric_expr_2x2(target_expr) + else: + expr = x2 + target_expr = x1 + + oneform_V = inner(expr, TestFunction(V)) * dx Q_target = make_quadrature_space(V) @@ -104,9 +116,8 @@ def test_cross_mesh_oneform_adjoint(V): # cofunc_V = assemble(interpolate(TestFunction(V_target), oneform_target)) # V^* -> V_target^* - cofunc_V_direct = assemble(inner(x1, TestFunction(V_target)) * dx) - + cofunc_V_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) assert np.allclose(cofunc_V.dat.data_ro, cofunc_V_direct.dat.data_ro) if __name__ == "__main__": - pytest.main([__file__ + "::test_cross_mesh_oneform_adjoint[RT_1]"]) \ No newline at end of file + pytest.main([__file__ + "::test_cross_mesh_oneform_adjoint[GLS_3]"]) \ No newline at end of file diff --git a/tsfc/driver.py b/tsfc/driver.py index d6528b779f..f83db1d944 100644 --- a/tsfc/driver.py +++ b/tsfc/driver.py @@ -245,20 +245,15 @@ def compile_expression_dual_evaluation(expression, ufl_element, *, complex_mode = is_complex(parameters["scalar_type"]) orig_coefficients = extract_coefficients(expression) - v, operand = expression.argument_slots() + if isinstance(expression, ufl.Interpolate): + v, operand = expression.argument_slots() + else: + operand = expression + v = ufl.FunctionSpace(extract_unique_domain(operand), ufl_element) # Map into reference space operand = apply_mapping(operand, ufl_element, domain) - if ufl_element.mapping() != "identity": - # Need to map dual argument for adjoint interpolation - ref_element = finat.ufl.WithMapping(ufl_element, "identity") - V = ufl.FunctionSpace(domain, ref_element) - if isinstance(v, ufl.Coargument): - v = ufl.Coargument(V.dual(), v.number()) - else: - v = ufl.Cofunction(V.dual()) - # Apply UFL preprocessing operand = ufl_utils.preprocess_expression(operand, complex_mode=complex_mode) operand = simplify_abs(operand, complex_mode) @@ -278,6 +273,9 @@ def compile_expression_dual_evaluation(expression, ufl_element, *, assert len(argument_multiindices) == len(arguments) # Replace coordinates (if any) unless otherwise specified by kwarg + if domain is None: + domain = extract_unique_domain(expression) + assert domain is not None builder._domain_integral_type_map = {domain: "cell"} builder._entity_ids = {domain: (0,)} From 277034046233508a384bfa3a1024f25abef6322a Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 21 Jan 2026 18:29:57 +0000 Subject: [PATCH 08/34] WIP --- firedrake/interpolation.py | 122 +++++++++++++++--- .../test_cross_mesh_non_lagrange.py | 79 +++++++++++- 2 files changed, 182 insertions(+), 19 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index b69fc22647..9576627502 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -28,7 +28,7 @@ from firedrake.utils import IntType, ScalarType, cached_property, known_pyop2_safe, tuplify from firedrake.pointeval_utils import runtime_quadrature_element from firedrake.tsfc_interface import extract_numbered_coefficients, _cachedir -from firedrake.ufl_expr import Argument, Coargument, action +from firedrake.ufl_expr import Argument, Coargument, TrialFunction, TestFunction, action from firedrake.mesh import MissingPointsBehaviour, VertexOnlyMeshTopology, MeshGeometry, MeshTopology, VertexOnlyMesh from firedrake.petsc import PETSc from firedrake.halo import _get_mtype @@ -430,7 +430,21 @@ def __init__(self, expr: Interpolate): if self.source_mesh.unique().geometric_dimension != self.target_mesh.unique().geometric_dimension: raise ValueError("Geometric dimensions of source and destination meshes must match.") - dest_element = self.target_space.ufl_element() + def _get_element(self, V: WithGeometry) -> FiniteElementBase: + """Return the element of the function space V. If V is tensor/vector valued, + return the base scalar element. + + Parameters + ---------- + V + A :class:`.WithGeometry` function space. + + Returns + ------- + FiniteElementBase + The base element of V. + """ + dest_element = V.ufl_element() if isinstance(dest_element, MixedElement): if isinstance(dest_element, VectorElement | TensorElement): # In this case all sub elements are equal @@ -440,16 +454,68 @@ def __init__(self, expr: Interpolate): "Can't yet cross-mesh interpolate onto function spaces made from VectorElements " "or TensorElements made from sub elements with value shape other than ()." ) - self.dest_element = base_element + return base_element else: raise NotImplementedError("Interpolation with MixedFunctionSpace requires MixedInterpolator.") else: # scalar fiat/finat element - self.dest_element = dest_element + return dest_element + + def _fs_type(self, V: WithGeometry) -> Callable[..., WithGeometry]: + """Returns a callable which returns a function space matching the type of V. + + Parameters + ---------- + V + A :class:`.WithGeometry` function space. + + Returns + ------- + Callable + A callable which returns a :class:`.WithGeometry` matching the type of V. + """ + # Get the correct type of function space + shape = V.ufl_function_space().value_shape + if len(shape) == 0: + return FunctionSpace + elif len(shape) == 1: + return partial(VectorFunctionSpace, dim=shape[0]) + else: + symmetry = V.ufl_element().symmetry() + return partial(TensorFunctionSpace, shape=shape, symmetry=symmetry) + + def _get_quadrature_space(self, V: WithGeometry) -> WithGeometry: + """Return a FunctionSpace whose element is a quadrature element with + points at the quadrature points of V's element. + + Used as an intermediate space for interpolating into a space which doesn't + have point evaluation dofs. + + Parameters + ---------- + V : WithGeometry + A :class:`.WithGeometry` function space to build the quadrature space for. + + Returns + ------- + WithGeometry + A :class:`.WithGeometry` function space with quadrature element. + """ + V_element = V.finat_element + _, points = V_element.dual_basis + weights = numpy.full(len(points.points), numpy.nan) # These can be any number since we never integrate + quad_scheme = QuadratureRule(points, weights, ref_el=V_element.cell) + element = FiniteElement("Quadrature", degree=V_element.degree, quad_scheme=quad_scheme) + return self._fs_type(V)(V.mesh(), element) + + def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpolate, Interpolate]: + """Return symbolic ``Interpolate`` expressions for point evaluation of the `target_space`s + dofs in the source mesh, and the corresponding input-ordering interpolation. - def _get_symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: - """Return the symbolic ``Interpolate`` expressions for point evaluation and - re-ordering into the input-ordering VertexOnlyMesh. + Parameters + ---------- + target_space + The :class:`.WithGeometry` function space which we are interpolating into. Returns ------- @@ -460,13 +526,19 @@ def _get_symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: Raises ------ DofNotDefinedError - If any DoFs in the target mesh cannot be defined in the source mesh. + If any of the target spaces dofs cannot be defined in the source mesh. + ValueError + If the target space does not have point-evaluation dofs. """ from firedrake.assemble import assemble + if target_space.ufl_element().mapping() != "identity": + raise ValueError(f"FunctionSpace {target_space} must have point-evaluation dofs.") + # Immerse coordinates of target space point evaluation dofs in src_mesh - target_space_vec = VectorFunctionSpace(self.target_mesh.unique(), self.dest_element) - f_dest_node_coords = assemble(interpolate(self.target_mesh.unique().coordinates, target_space_vec)) - dest_node_coords = f_dest_node_coords.dat.data_ro.reshape(-1, self.target_mesh.unique().geometric_dimension) + target_mesh = target_space.mesh().unique() + target_space_vec = VectorFunctionSpace(target_mesh, self._get_element(target_space)) + f_dest_node_coords = assemble(interpolate(target_mesh.coordinates, target_space_vec)) + dest_node_coords = f_dest_node_coords.dat.data_ro.reshape(-1, target_mesh.geometric_dimension) try: vom = VertexOnlyMesh( self.source_mesh.unique(), @@ -492,6 +564,7 @@ def _get_symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: fs_type = partial(TensorFunctionSpace, shape=shape, symmetry=symmetry) # Get expression for point evaluation at the dest_node_coords + fs_type = self._fs_type(target_space) P0DG_vom = fs_type(vom, "DG", 0) point_eval = interpolate(self.operand, P0DG_vom) @@ -510,9 +583,16 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) # self.ufl_interpolate.function_space() is None in the 0-form case V_dest = self.ufl_interpolate.function_space() or self.target_space - f = tensor or Function(V_dest) - point_eval, point_eval_input_ordering = self._get_symbolic_expressions() + # Interpolate into intermediate quadrature space for non-identity mapped elements + if into_quadrature_space := V_dest.ufl_element().mapping() != "identity": + Q_dest = self._get_quadrature_space(V_dest) + else: + Q_dest = None + target_space = Q_dest or V_dest + f = tensor or Function(target_space) + + point_eval, point_eval_input_ordering = self._get_symbolic_expressions(target_space) P0DG_vom_input_ordering = point_eval_input_ordering.argument_slots()[0].function_space().dual() if self.rank == 2: @@ -521,12 +601,22 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) # `self.point_eval_interpolate` and the permutation # given by `self.to_input_ordering_interpolate`. if self.ufl_interpolate.is_adjoint: - symbolic = action(point_eval, point_eval_input_ordering) + interp_expr = action(point_eval, point_eval_input_ordering) else: - symbolic = action(point_eval_input_ordering, point_eval) + interp_expr = action(point_eval_input_ordering, point_eval) def callable() -> PETSc.Mat: - return assemble(symbolic, mat_type=mat_type).petscmat + res = assemble(interp_expr, mat_type=mat_type).petscmat + if into_quadrature_space: + if self.ufl_interpolate.is_adjoint: + I = AssembledMatrix((Argument(Q_dest, 0), Argument(V_dest.dual(), 1)), None, res) + return assemble(action(interpolate(TestFunction(V_dest), Q_dest))).petscmat + else: + I = AssembledMatrix((Argument(Q_dest.dual(), 0), Argument(V_dest, 1)), None, res) # V_dest x Q_dest^* -> R + return assemble(action(interpolate(TrialFunction(Q_dest), V_dest), I)).petscmat # Q_dest x V_dest^* -> R + else: + return res + elif self.ufl_interpolate.is_adjoint: assert self.rank == 1 # f_src is a cofunction on V_dest.dual diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 6651eada43..3d0e1b0e35 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -70,11 +70,11 @@ def test_cross_mesh_oneform(V): Q_target = make_quadrature_space(V) # Interp V_source -> Q - I1 = interpolate(f_source, Q_target) + I1 = interpolate(f_source, Q_target) # CrossMeshInterpolator f_quadrature = assemble(I1) # Interp Q -> V_target - I2 = interpolate(f_quadrature, V) + I2 = interpolate(f_quadrature, V) # SameMeshInterpolator f_target = assemble(I2) f_direct = Function(V).interpolate(expr2) @@ -82,6 +82,42 @@ def test_cross_mesh_oneform(V): assert np.allclose(f_target.dat.data_ro, f_direct.dat.data_ro) +def test_cross_mesh_twoform(V): + mesh1 = UnitSquareMesh(5, 5) + mesh2 = V.mesh() + x, y = SpatialCoordinate(mesh1) + x1, y1 = SpatialCoordinate(mesh2) + + shape = V.ufl_function_space().value_shape + if len(shape) == 1: + fs_type = partial(VectorFunctionSpace, dim=shape[0]) + expr1 = as_vector([x, y]) + expr2 = as_vector([x1, y1]) + elif len(shape) == 2: + fs_type = partial(TensorFunctionSpace, shape=shape) + expr1 = as_tensor([[x, x*y], [x*y, y]]) + expr2 = as_tensor([[x1, x1*y1], [x1*y1, y1]]) + else: + raise ValueError("Unsupported target space shape") + + V_source = fs_type(mesh1, "CG", 2) + Q_target = make_quadrature_space(V) + + I1 = interpolate(TrialFunction(V_source), Q_target) # V_source x Q_target^* -> R + I2 = interpolate(TrialFunction(Q_target), V) # Q_target x V^* -> R + + I = assemble(action(I2, I1)) # V_source x V^* -> R + + I_direct = assemble(interpolate(TrialFunction(V_source), V)) # V_source x V^* -> R + + f_source = Function(V_source).interpolate(expr1) + f_direct = Function(V).interpolate(expr2) + + f_interpolated = assemble(action(I_direct, f_source)) + assert np.allclose(f_interpolated.dat.data_ro, f_direct.dat.data_ro) + + + def test_cross_mesh_oneform_adjoint(V): # Can already do Lagrange -> RT adjoint # V^* -> Q^* -> V_target^* @@ -119,5 +155,42 @@ def test_cross_mesh_oneform_adjoint(V): cofunc_V_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) assert np.allclose(cofunc_V.dat.data_ro, cofunc_V_direct.dat.data_ro) + +def test_cross_mesh_twoform_adjoint(V): + # V^* -> Q^* -> V_target^* + mesh1 = UnitSquareMesh(2, 2) + x1 = SpatialCoordinate(mesh1) + V_target = fs_shape(V)(mesh1, "CG", 1) + mesh2 = V.mesh() + x2 = SpatialCoordinate(mesh2) + + if len(V.value_shape) > 1: + expr = outer(x2, x2) + target_expr = outer(x1, x1) + if V.ufl_element().mapping() == "covariant contravariant Piola": + expr = deviatoric_expr_2x2(expr) + target_expr = deviatoric_expr_2x2(target_expr) + else: + expr = x2 + target_expr = x1 + + oneform_V = inner(expr, TestFunction(V)) * dx + + Q = make_quadrature_space(V) + + I1 = interpolate(TestFunction(Q), V) # V^* x Q -> R + I2 = interpolate(TestFunction(V_target), Q) # Q^* x V_target -> R + I = assemble(action(I2, I1)) # V^* x V_target -> R + + # I_direct = assemble(interpolate(TestFunction(V_target), V)) # V^* x V_target -> R + # assert I_direct.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) + + cofunc_V = assemble(action(I, oneform_V)) + cofunc_V_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) + + assert np.allclose(cofunc_V.dat.data_ro, cofunc_V_direct.dat.data_ro) + + + if __name__ == "__main__": - pytest.main([__file__ + "::test_cross_mesh_oneform_adjoint[GLS_3]"]) \ No newline at end of file + pytest.main([__file__ + "::test_cross_mesh_twoform[RT_2]"]) \ No newline at end of file From cdf2a8c901fb8b3e9bcb490b5d1518dd618fc5ae Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 21 Jan 2026 18:53:50 +0000 Subject: [PATCH 09/34] fix target_space --- firedrake/interpolation.py | 15 ++++++++------- .../regression/test_cross_mesh_non_lagrange.py | 6 +++--- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 9576627502..abd6b75792 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -583,14 +583,15 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) # self.ufl_interpolate.function_space() is None in the 0-form case V_dest = self.ufl_interpolate.function_space() or self.target_space + f = tensor or Function(V_dest) # Interpolate into intermediate quadrature space for non-identity mapped elements - if into_quadrature_space := V_dest.ufl_element().mapping() != "identity": + if into_quadrature_space := self.target_space.ufl_element().mapping() != "identity": Q_dest = self._get_quadrature_space(V_dest) else: Q_dest = None - target_space = Q_dest or V_dest - f = tensor or Function(target_space) + target_space = Q_dest or self.target_space + point_eval, point_eval_input_ordering = self._get_symbolic_expressions(target_space) P0DG_vom_input_ordering = point_eval_input_ordering.argument_slots()[0].function_space().dual() @@ -609,11 +610,11 @@ def callable() -> PETSc.Mat: res = assemble(interp_expr, mat_type=mat_type).petscmat if into_quadrature_space: if self.ufl_interpolate.is_adjoint: - I = AssembledMatrix((Argument(Q_dest, 0), Argument(V_dest.dual(), 1)), None, res) - return assemble(action(interpolate(TestFunction(V_dest), Q_dest))).petscmat + I = AssembledMatrix((Argument(Q_dest, 0), Argument(target_space.dual(), 1)), None, res) + return assemble(action(interpolate(TestFunction(target_space), Q_dest), I)).petscmat else: - I = AssembledMatrix((Argument(Q_dest.dual(), 0), Argument(V_dest, 1)), None, res) # V_dest x Q_dest^* -> R - return assemble(action(interpolate(TrialFunction(Q_dest), V_dest), I)).petscmat # Q_dest x V_dest^* -> R + I = AssembledMatrix((Argument(Q_dest.dual(), 0), Argument(target_space, 1)), None, res) # target_space x Q_dest^* -> R + return assemble(action(interpolate(TrialFunction(Q_dest), target_space), I)).petscmat # Q_dest x target_space^* -> R else: return res diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 3d0e1b0e35..114470d83a 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -182,10 +182,10 @@ def test_cross_mesh_twoform_adjoint(V): I2 = interpolate(TestFunction(V_target), Q) # Q^* x V_target -> R I = assemble(action(I2, I1)) # V^* x V_target -> R - # I_direct = assemble(interpolate(TestFunction(V_target), V)) # V^* x V_target -> R - # assert I_direct.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) + I_direct = assemble(interpolate(TestFunction(V_target), V)) # V^* x V_target -> R + assert I_direct.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) - cofunc_V = assemble(action(I, oneform_V)) + cofunc_V = assemble(action(I_direct, oneform_V)) cofunc_V_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) assert np.allclose(cofunc_V.dat.data_ro, cofunc_V_direct.dat.data_ro) From 4de6baf8a940b8bcc518fb322850027c895695f7 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 21 Jan 2026 19:14:13 +0000 Subject: [PATCH 10/34] two-form assembly working --- firedrake/interpolation.py | 8 ++++---- .../firedrake/regression/test_cross_mesh_non_lagrange.py | 3 ++- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index abd6b75792..cabe230755 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -587,7 +587,7 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) # Interpolate into intermediate quadrature space for non-identity mapped elements if into_quadrature_space := self.target_space.ufl_element().mapping() != "identity": - Q_dest = self._get_quadrature_space(V_dest) + Q_dest = self._get_quadrature_space(self.target_space) else: Q_dest = None target_space = Q_dest or self.target_space @@ -611,10 +611,10 @@ def callable() -> PETSc.Mat: if into_quadrature_space: if self.ufl_interpolate.is_adjoint: I = AssembledMatrix((Argument(Q_dest, 0), Argument(target_space.dual(), 1)), None, res) - return assemble(action(interpolate(TestFunction(target_space), Q_dest), I)).petscmat + return assemble(action(I, interpolate(TestFunction(Q_dest), self.target_space))).petscmat else: - I = AssembledMatrix((Argument(Q_dest.dual(), 0), Argument(target_space, 1)), None, res) # target_space x Q_dest^* -> R - return assemble(action(interpolate(TrialFunction(Q_dest), target_space), I)).petscmat # Q_dest x target_space^* -> R + I = AssembledMatrix((Argument(Q_dest.dual(), 0), Argument(target_space, 1)), None, res) + return assemble(action(interpolate(TrialFunction(Q_dest), self.target_space), I)).petscmat else: return res diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 114470d83a..4e9cd6f122 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -181,6 +181,7 @@ def test_cross_mesh_twoform_adjoint(V): I1 = interpolate(TestFunction(Q), V) # V^* x Q -> R I2 = interpolate(TestFunction(V_target), Q) # Q^* x V_target -> R I = assemble(action(I2, I1)) # V^* x V_target -> R + assert I.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) I_direct = assemble(interpolate(TestFunction(V_target), V)) # V^* x V_target -> R assert I_direct.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) @@ -193,4 +194,4 @@ def test_cross_mesh_twoform_adjoint(V): if __name__ == "__main__": - pytest.main([__file__ + "::test_cross_mesh_twoform[RT_2]"]) \ No newline at end of file + pytest.main([__file__ + "::test_cross_mesh_twoform_adjoint[RT_2]"]) \ No newline at end of file From c8f6e08ac968b31f7dffc461a35fab697a6ead67 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 26 Jan 2026 16:02:09 +0000 Subject: [PATCH 11/34] forward and adjoint one-form --- firedrake/interpolation.py | 50 +++++++------- .../test_cross_mesh_non_lagrange.py | 66 ++++++------------- .../regression/test_interpolate_cross_mesh.py | 1 - 3 files changed, 44 insertions(+), 73 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index cabe230755..268c9efda2 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -475,7 +475,7 @@ def _fs_type(self, V: WithGeometry) -> Callable[..., WithGeometry]: A callable which returns a :class:`.WithGeometry` matching the type of V. """ # Get the correct type of function space - shape = V.ufl_function_space().value_shape + shape = V.value_shape if len(shape) == 0: return FunctionSpace elif len(shape) == 1: @@ -485,7 +485,7 @@ def _fs_type(self, V: WithGeometry) -> Callable[..., WithGeometry]: return partial(TensorFunctionSpace, shape=shape, symmetry=symmetry) def _get_quadrature_space(self, V: WithGeometry) -> WithGeometry: - """Return a FunctionSpace whose element is a quadrature element with + """Return a FunctionSpace whose element is a quadrature element with points at the quadrature points of V's element. Used as an intermediate space for interpolating into a space which doesn't @@ -500,13 +500,12 @@ def _get_quadrature_space(self, V: WithGeometry) -> WithGeometry: ------- WithGeometry A :class:`.WithGeometry` function space with quadrature element. - """ + """ V_element = V.finat_element _, points = V_element.dual_basis weights = numpy.full(len(points.points), numpy.nan) # These can be any number since we never integrate quad_scheme = QuadratureRule(points, weights, ref_el=V_element.cell) - element = FiniteElement("Quadrature", degree=V_element.degree, quad_scheme=quad_scheme) - return self._fs_type(V)(V.mesh(), element) + return self._fs_type(V)(V.mesh(), "Quadrature", degree=V_element.degree, quad_scheme=quad_scheme) def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpolate, Interpolate]: """Return symbolic ``Interpolate`` expressions for point evaluation of the `target_space`s @@ -553,16 +552,6 @@ def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpo "This may be because the target mesh covers a larger domain than the " "source mesh. To disable this error, set allow_missing_dofs=True.") - # Get the correct type of function space - shape = self.target_space.ufl_function_space().value_shape - if len(shape) == 0: - fs_type = FunctionSpace - elif len(shape) == 1: - fs_type = partial(VectorFunctionSpace, dim=shape[0]) - else: - symmetry = self.target_space.ufl_element().symmetry() - fs_type = partial(TensorFunctionSpace, shape=shape, symmetry=symmetry) - # Get expression for point evaluation at the dest_node_coords fs_type = self._fs_type(target_space) P0DG_vom = fs_type(vom, "DG", 0) @@ -581,17 +570,15 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) raise NotImplementedError("bcs not implemented for cross-mesh interpolation.") mat_type = mat_type or "aij" - # self.ufl_interpolate.function_space() is None in the 0-form case - V_dest = self.ufl_interpolate.function_space() or self.target_space - f = tensor or Function(V_dest) - # Interpolate into intermediate quadrature space for non-identity mapped elements if into_quadrature_space := self.target_space.ufl_element().mapping() != "identity": Q_dest = self._get_quadrature_space(self.target_space) else: Q_dest = None target_space = Q_dest or self.target_space - + + # self.ufl_interpolate.function_space() is None in the 0-form case + f = tensor or Function(target_space.dual() if self.ufl_interpolate.is_adjoint else target_space) point_eval, point_eval_input_ordering = self._get_symbolic_expressions(target_space) P0DG_vom_input_ordering = point_eval_input_ordering.argument_slots()[0].function_space().dual() @@ -620,8 +607,12 @@ def callable() -> PETSc.Mat: elif self.ufl_interpolate.is_adjoint: assert self.rank == 1 - # f_src is a cofunction on V_dest.dual - cofunc = self.dual_arg + + if into_quadrature_space: + cofunc = assemble(interpolate(TestFunction(target_space), self.dual_arg)) + else: + cofunc = self.dual_arg + assert isinstance(cofunc, Cofunction) # Our first adjoint operation is to assign the dat values to a @@ -636,8 +627,9 @@ def callable() -> PETSc.Mat: # and all points from the input ordering VOM are in the original. def callable() -> Cofunction: f_src_at_src_node_coords = assemble(action(point_eval_input_ordering, f_input_ordering)) - assemble(action(point_eval, f_src_at_src_node_coords), tensor=f) - return f + f_target = Cofunction(point_eval.function_space()) if into_quadrature_space else f + assemble(action(point_eval, f_src_at_src_node_coords), tensor=f_target) + return f_target else: assert self.rank in {0, 1} # We create the input-ordering Function before interpolating so we can @@ -661,12 +653,18 @@ def callable() -> Function | Number: else: f.dat.data_wo[:] = f_point_eval_input_ordering.dat.data_ro[:] + if into_quadrature_space: + f_target = Function(self.target_space) + assemble(interpolate(f, self.target_space), tensor=f_target) + else: + f_target = f + if self.rank == 0: # We take the action of the dual_arg on the interpolated function assert isinstance(self.dual_arg, Cofunction) - return assemble(action(self.dual_arg, f)) + return assemble(action(self.dual_arg, f_target)) else: - return f + return f_target return callable @property diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 4e9cd6f122..837d3ce1eb 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -17,9 +17,10 @@ def fs_shape(V): raise ValueError("Invalid function space shape") return fs_type + def make_quadrature_space(V): """Builds a quadrature space on the target mesh. - + This space has point evaluation dofs at the quadrature PointSet of the target space's element @@ -33,17 +34,15 @@ def make_quadrature_space(V): @pytest.fixture(params=[("RT", 2), ("RT", 3), ("RT", 4), ("BDM", 1), ("BDM", 2), ("BDM", 3), - ("BDFM", 2), ("HHJ", 2),("N1curl", 2), ("N1curl", 3), ("N1curl", 4), + ("BDFM", 2), ("HHJ", 2), ("N1curl", 2), ("N1curl", 3), ("N1curl", 4), ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 2), ("GLS", 3), ("GLS", 4), ("GLS2", 1), ("GLS2", 2), ("GLS2", 3)], - ids=lambda x: f"{x[0]}_{x[1]}") + ids=lambda x: f"{x[0]}_{x[1]}") def V(request): element, degree = request.param mesh = UnitSquareMesh(8, 8) return FunctionSpace(mesh, element, degree) -# V_source -> Q -> V_target -# V_target^* -> Q^* -> V_source^* def test_cross_mesh_oneform(V): mesh1 = UnitSquareMesh(5, 5) @@ -62,21 +61,11 @@ def test_cross_mesh_oneform(V): expr2 = as_tensor([[x1, x1*y1], [x1*y1, y1]]) else: raise ValueError("Unsupported target space shape") - + V_source = fs_type(mesh1, "CG", 2) f_source = Function(V_source).interpolate(expr1) - # Make intermediate Quadrature space on target mesh - Q_target = make_quadrature_space(V) - - # Interp V_source -> Q - I1 = interpolate(f_source, Q_target) # CrossMeshInterpolator - f_quadrature = assemble(I1) - - # Interp Q -> V_target - I2 = interpolate(f_quadrature, V) # SameMeshInterpolator - f_target = assemble(I2) - + f_target = assemble(interpolate(f_source, V)) f_direct = Function(V).interpolate(expr2) assert np.allclose(f_target.dat.data_ro, f_direct.dat.data_ro) @@ -99,25 +88,18 @@ def test_cross_mesh_twoform(V): expr2 = as_tensor([[x1, x1*y1], [x1*y1, y1]]) else: raise ValueError("Unsupported target space shape") - - V_source = fs_type(mesh1, "CG", 2) - Q_target = make_quadrature_space(V) - I1 = interpolate(TrialFunction(V_source), Q_target) # V_source x Q_target^* -> R - I2 = interpolate(TrialFunction(Q_target), V) # Q_target x V^* -> R - - I = assemble(action(I2, I1)) # V_source x V^* -> R + V_source = fs_type(mesh1, "CG", 2) - I_direct = assemble(interpolate(TrialFunction(V_source), V)) # V_source x V^* -> R + I = assemble(interpolate(TrialFunction(V_source), V)) # V_source x V^* -> R f_source = Function(V_source).interpolate(expr1) f_direct = Function(V).interpolate(expr2) - f_interpolated = assemble(action(I_direct, f_source)) + f_interpolated = assemble(action(I, f_source)) assert np.allclose(f_interpolated.dat.data_ro, f_direct.dat.data_ro) - def test_cross_mesh_oneform_adjoint(V): # Can already do Lagrange -> RT adjoint # V^* -> Q^* -> V_target^* @@ -140,17 +122,17 @@ def test_cross_mesh_oneform_adjoint(V): oneform_V = inner(expr, TestFunction(V)) * dx - Q_target = make_quadrature_space(V) + # Q_target = make_quadrature_space(V) - # Interp V^* -> Q^* - I1_adj = interpolate(TestFunction(Q_target), oneform_V) - cofunc_Q = assemble(I1_adj) + # # Interp V^* -> Q^* + # I1_adj = interpolate(TestFunction(Q_target), oneform_V) # SameMesh + # cofunc_Q = assemble(I1_adj) - # Interp Q^* -> V_target^* - I2_adj = interpolate(TestFunction(V_target), cofunc_Q) - cofunc_V = assemble(I2_adj) + # # Interp Q^* -> V_target^* + # I2_adj = interpolate(TestFunction(V_target), cofunc_Q) # CrossMesh + # cofunc_V = assemble(I2_adj) - # cofunc_V = assemble(interpolate(TestFunction(V_target), oneform_target)) # V^* -> V_target^* + cofunc_V = assemble(interpolate(TestFunction(V_target), oneform_V)) # V^* -> V_target^* cofunc_V_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) assert np.allclose(cofunc_V.dat.data_ro, cofunc_V_direct.dat.data_ro) @@ -176,22 +158,14 @@ def test_cross_mesh_twoform_adjoint(V): oneform_V = inner(expr, TestFunction(V)) * dx - Q = make_quadrature_space(V) - - I1 = interpolate(TestFunction(Q), V) # V^* x Q -> R - I2 = interpolate(TestFunction(V_target), Q) # Q^* x V_target -> R - I = assemble(action(I2, I1)) # V^* x V_target -> R + I = assemble(interpolate(TestFunction(V_target), V)) # V^* x V_target -> R assert I.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) - I_direct = assemble(interpolate(TestFunction(V_target), V)) # V^* x V_target -> R - assert I_direct.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) - - cofunc_V = assemble(action(I_direct, oneform_V)) + cofunc_V = assemble(action(I, oneform_V)) cofunc_V_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) assert np.allclose(cofunc_V.dat.data_ro, cofunc_V_direct.dat.data_ro) - if __name__ == "__main__": - pytest.main([__file__ + "::test_cross_mesh_twoform_adjoint[RT_2]"]) \ No newline at end of file + pytest.main([__file__ + "::test_cross_mesh_oneform_adjoint[RT_2]"]) diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index c57a8e025e..c256f15f69 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -6,7 +6,6 @@ import numpy as np import pytest from ufl import product -import subprocess def allgather(comm, coords): From 97ce4ea8ed43b3ddb60270571d0252cfb0cd9cb4 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 26 Jan 2026 16:19:35 +0000 Subject: [PATCH 12/34] fixes --- firedrake/interpolation.py | 35 ++++++++++--------- .../test_cross_mesh_non_lagrange.py | 21 +++++------ 2 files changed, 29 insertions(+), 27 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 268c9efda2..4bc468f628 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -578,7 +578,7 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) target_space = Q_dest or self.target_space # self.ufl_interpolate.function_space() is None in the 0-form case - f = tensor or Function(target_space.dual() if self.ufl_interpolate.is_adjoint else target_space) + f = tensor or Function(self.ufl_interpolate.function_space() or self.target_space) point_eval, point_eval_input_ordering = self._get_symbolic_expressions(target_space) P0DG_vom_input_ordering = point_eval_input_ordering.argument_slots()[0].function_space().dual() @@ -608,26 +608,27 @@ def callable() -> PETSc.Mat: elif self.ufl_interpolate.is_adjoint: assert self.rank == 1 - if into_quadrature_space: - cofunc = assemble(interpolate(TestFunction(target_space), self.dual_arg)) - else: - cofunc = self.dual_arg + def callable() -> Cofunction: + if into_quadrature_space: + cofunc = assemble(interpolate(TestFunction(target_space), self.dual_arg)) + f_target = Cofunction(point_eval.function_space()) + else: + cofunc = self.dual_arg + f_target = f - assert isinstance(cofunc, Cofunction) + assert isinstance(cofunc, Cofunction) - # Our first adjoint operation is to assign the dat values to a - # P0DG cofunction on our input ordering VOM. - f_input_ordering = Cofunction(P0DG_vom_input_ordering.dual()) - f_input_ordering.dat.data_wo[:] = cofunc.dat.data_ro[:] + # Our first adjoint operation is to assign the dat values to a + # P0DG cofunction on our input ordering VOM. + f_input_ordering = Cofunction(P0DG_vom_input_ordering.dual()) + f_input_ordering.dat.data_wo[:] = cofunc.dat.data_ro[:] - # The rest of the adjoint interpolation is the composition - # of the adjoint interpolators in the reverse direction. - # We don't worry about skipping over missing points here - # because we're going from the input ordering VOM to the original VOM - # and all points from the input ordering VOM are in the original. - def callable() -> Cofunction: + # The rest of the adjoint interpolation is the composition + # of the adjoint interpolators in the reverse direction. + # We don't worry about skipping over missing points here + # because we're going from the input ordering VOM to the original VOM + # and all points from the input ordering VOM are in the original. f_src_at_src_node_coords = assemble(action(point_eval_input_ordering, f_input_ordering)) - f_target = Cofunction(point_eval.function_space()) if into_quadrature_space else f assemble(action(point_eval, f_src_at_src_node_coords), tensor=f_target) return f_target else: diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 837d3ce1eb..d5155a7b40 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -122,20 +122,21 @@ def test_cross_mesh_oneform_adjoint(V): oneform_V = inner(expr, TestFunction(V)) * dx - # Q_target = make_quadrature_space(V) + Q_target = make_quadrature_space(V) - # # Interp V^* -> Q^* - # I1_adj = interpolate(TestFunction(Q_target), oneform_V) # SameMesh - # cofunc_Q = assemble(I1_adj) + # Interp V^* -> Q^* + I1_adj = interpolate(TestFunction(Q_target), oneform_V) # SameMesh + cofunc_Q = assemble(I1_adj) - # # Interp Q^* -> V_target^* - # I2_adj = interpolate(TestFunction(V_target), cofunc_Q) # CrossMesh - # cofunc_V = assemble(I2_adj) + # Interp Q^* -> V_target^* + I2_adj = interpolate(TestFunction(V_target), cofunc_Q) # CrossMesh + cofunc_Vtarget_manual = assemble(I2_adj) - cofunc_V = assemble(interpolate(TestFunction(V_target), oneform_V)) # V^* -> V_target^* + cofunc_Vtarget = assemble(interpolate(TestFunction(V_target), oneform_V)) # V^* -> V_target^* + assert np.allclose(cofunc_Vtarget_manual.dat.data_ro, cofunc_Vtarget.dat.data_ro) - cofunc_V_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) - assert np.allclose(cofunc_V.dat.data_ro, cofunc_V_direct.dat.data_ro) + cofunc_Vtarget_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) + assert np.allclose(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) def test_cross_mesh_twoform_adjoint(V): From b2f3a876ce970d6f49ddbfa64379220d405d1d38 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 26 Jan 2026 23:07:28 +0000 Subject: [PATCH 13/34] working remove --- firedrake/interpolation.py | 8 +- .../test_cross_mesh_non_lagrange.py | 152 ++++++++---------- .../regression/test_interpolate_cross_mesh.py | 14 +- 3 files changed, 78 insertions(+), 96 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 4bc468f628..c853ac326f 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -577,8 +577,12 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) Q_dest = None target_space = Q_dest or self.target_space - # self.ufl_interpolate.function_space() is None in the 0-form case - f = tensor or Function(self.ufl_interpolate.function_space() or self.target_space) + if into_quadrature_space and not self.ufl_interpolate.is_adjoint: + f = Function(target_space) + else: + # self.ufl_interpolate.function_space() is None in the 0-form case + f = Function(self.ufl_interpolate.function_space() or target_space) + # f = tensor or Function(target_space.dual() if self.ufl_interpolate.is_adjoint else target_space) point_eval, point_eval_input_ordering = self._get_symbolic_expressions(target_space) P0DG_vom_input_ordering = point_eval_input_ordering.argument_slots()[0].function_space().dual() diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index d5155a7b40..4796cb2452 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -4,8 +4,6 @@ from finat.quadrature import QuadratureRule from functools import partial -from ufl.compound_expressions import deviatoric_expr_2x2 - def fs_shape(V): shape = V.ufl_function_space().value_shape @@ -44,34 +42,8 @@ def V(request): return FunctionSpace(mesh, element, degree) -def test_cross_mesh_oneform(V): - mesh1 = UnitSquareMesh(5, 5) - mesh2 = V.mesh() - x, y = SpatialCoordinate(mesh1) - x1, y1 = SpatialCoordinate(mesh2) - - shape = V.ufl_function_space().value_shape - if len(shape) == 1: - fs_type = partial(VectorFunctionSpace, dim=shape[0]) - expr1 = as_vector([x, y]) - expr2 = as_vector([x1, y1]) - elif len(shape) == 2: - fs_type = partial(TensorFunctionSpace, shape=shape) - expr1 = as_tensor([[x, x*y], [x*y, y]]) - expr2 = as_tensor([[x1, x1*y1], [x1*y1, y1]]) - else: - raise ValueError("Unsupported target space shape") - - V_source = fs_type(mesh1, "CG", 2) - f_source = Function(V_source).interpolate(expr1) - - f_target = assemble(interpolate(f_source, V)) - f_direct = Function(V).interpolate(expr2) - - assert np.allclose(f_target.dat.data_ro, f_direct.dat.data_ro) - - -def test_cross_mesh_twoform(V): +@pytest.mark.parametrize("rank", [1, 2]) +def test_cross_mesh(V, rank): mesh1 = UnitSquareMesh(5, 5) mesh2 = V.mesh() x, y = SpatialCoordinate(mesh1) @@ -90,17 +62,42 @@ def test_cross_mesh_twoform(V): raise ValueError("Unsupported target space shape") V_source = fs_type(mesh1, "CG", 2) - - I = assemble(interpolate(TrialFunction(V_source), V)) # V_source x V^* -> R - f_source = Function(V_source).interpolate(expr1) f_direct = Function(V).interpolate(expr2) - f_interpolated = assemble(action(I, f_source)) - assert np.allclose(f_interpolated.dat.data_ro, f_direct.dat.data_ro) - - -def test_cross_mesh_oneform_adjoint(V): + Q = make_quadrature_space(V) + + if rank == 2: + # Assemble the operator + I1 = interpolate(TrialFunction(V_source), Q) # V_source x Q_target^* -> R + I2 = interpolate(TrialFunction(Q), V) # Q_target x V^* -> R + I_manual = assemble(action(I2, I1)) # V_source x V^* -> R + assert I_manual.arguments() == (TestFunction(V.dual()), TrialFunction(V_source)) + # Direct assembly + I_direct = assemble(interpolate(TrialFunction(V_source), V)) # V_source + assert I_direct.arguments() == (TestFunction(V.dual()), TrialFunction(V_source)) + + f_interpolated_manual = assemble(action(I_manual, f_source)) + assert np.allclose(f_interpolated_manual.dat.data_ro, f_direct.dat.data_ro) + f_interpolated_direct = assemble(action(I_direct, f_source)) + assert np.allclose(f_interpolated_direct.dat.data_ro, f_direct.dat.data_ro) + elif rank == 1: + # Interp V_source -> Q + I1 = interpolate(f_source, Q) # SameMesh + f_quadrature = assemble(I1) + # Interp Q -> V + I2 = interpolate(f_quadrature, V) # CrossMesh + f_interpolated_manual = assemble(I2) + assert f_interpolated_manual.function_space() == V + assert np.allclose(f_interpolated_manual.dat.data_ro, f_direct.dat.data_ro) + + f_interpolated_direct = assemble(interpolate(f_source, V)) + assert f_interpolated_direct.function_space() == V + assert np.allclose(f_interpolated_direct.dat.data_ro, f_direct.dat.data_ro) + + +@pytest.mark.parametrize("rank", [1, 2]) +def test_cross_mesh_adjoint(V, rank): # Can already do Lagrange -> RT adjoint # V^* -> Q^* -> V_target^* mesh1 = UnitSquareMesh(2, 2) @@ -114,59 +111,40 @@ def test_cross_mesh_oneform_adjoint(V): expr = outer(x2, x2) target_expr = outer(x1, x1) if V.ufl_element().mapping() == "covariant contravariant Piola": - expr = deviatoric_expr_2x2(expr) - target_expr = deviatoric_expr_2x2(target_expr) + expr = dev(expr) + target_expr = dev(target_expr) else: expr = x2 target_expr = x1 oneform_V = inner(expr, TestFunction(V)) * dx - - Q_target = make_quadrature_space(V) - - # Interp V^* -> Q^* - I1_adj = interpolate(TestFunction(Q_target), oneform_V) # SameMesh - cofunc_Q = assemble(I1_adj) - - # Interp Q^* -> V_target^* - I2_adj = interpolate(TestFunction(V_target), cofunc_Q) # CrossMesh - cofunc_Vtarget_manual = assemble(I2_adj) - - cofunc_Vtarget = assemble(interpolate(TestFunction(V_target), oneform_V)) # V^* -> V_target^* - assert np.allclose(cofunc_Vtarget_manual.dat.data_ro, cofunc_Vtarget.dat.data_ro) - cofunc_Vtarget_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) - assert np.allclose(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) - - -def test_cross_mesh_twoform_adjoint(V): - # V^* -> Q^* -> V_target^* - mesh1 = UnitSquareMesh(2, 2) - x1 = SpatialCoordinate(mesh1) - V_target = fs_shape(V)(mesh1, "CG", 1) - mesh2 = V.mesh() - x2 = SpatialCoordinate(mesh2) - - if len(V.value_shape) > 1: - expr = outer(x2, x2) - target_expr = outer(x1, x1) - if V.ufl_element().mapping() == "covariant contravariant Piola": - expr = deviatoric_expr_2x2(expr) - target_expr = deviatoric_expr_2x2(target_expr) - else: - expr = x2 - target_expr = x1 - - oneform_V = inner(expr, TestFunction(V)) * dx - - I = assemble(interpolate(TestFunction(V_target), V)) # V^* x V_target -> R - assert I.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) - - cofunc_V = assemble(action(I, oneform_V)) - cofunc_V_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) - - assert np.allclose(cofunc_V.dat.data_ro, cofunc_V_direct.dat.data_ro) - -if __name__ == "__main__": - pytest.main([__file__ + "::test_cross_mesh_oneform_adjoint[RT_2]"]) + Q = make_quadrature_space(V) + + if rank == 2: + # Assemble the operator + I1 = interpolate(TestFunction(Q), V) # V^* x Q -> R + I2 = interpolate(TestFunction(V_target), Q) # Q^* x V_target -> R + I_manual = assemble(action(I2, I1)) # V^* x V_target -> R + assert I_manual.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) + # Direct assembly + I_direct = assemble(interpolate(TestFunction(V_target), V)) # V^* x V_target -> R + assert I_direct.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) + + cofunc_Vtarget_manual = assemble(action(I_manual, oneform_V)) + assert np.allclose(cofunc_Vtarget_manual.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) + cofunc_Vtarget = assemble(action(I_direct, oneform_V)) + assert np.allclose(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) + elif rank == 1: + # Interp V^* -> Q^* + I1_adj = interpolate(TestFunction(Q), oneform_V) # SameMesh + cofunc_Q = assemble(I1_adj) + + # Interp Q^* -> V_target^* + I2_adj = interpolate(TestFunction(V_target), cofunc_Q) # CrossMesh + cofunc_Vtarget_manual = assemble(I2_adj) + assert np.allclose(cofunc_Vtarget_manual.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) + + cofunc_Vtarget = assemble(interpolate(TestFunction(V_target), oneform_V)) # V^* -> V_target^* + assert np.allclose(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index c256f15f69..b721fb40d2 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -16,9 +16,9 @@ def allgather(comm, coords): return coords -def unitsquaresetup(): +def unitsquaresetup(dest_quad=True): m_src = UnitSquareMesh(2, 3) - m_dest = UnitSquareMesh(3, 5, quadrilateral=True) + m_dest = UnitSquareMesh(3, 5, quadrilateral=dest_quad) coords = np.array( [[0.56, 0.6], [0.1, 0.9], [0.9, 0.1], [0.9, 0.9], [0.726, 0.6584]] ) # fairly arbitrary @@ -39,7 +39,7 @@ def make_high_order(m_low_order, degree): @pytest.fixture( params=[ - "unitsquare", + "unitsquare_RT_N1curl_destination", "circlemanifold", "circlemanifold_to_high_order", "unitsquare_from_high_order", @@ -48,7 +48,7 @@ def make_high_order(m_low_order, degree): "unitsquare_vfs", "unitsquare_tfs", "unitsquare_N1curl_source", - "unitsquare_SminusDiv_destination", + "unitsquare_RT_N1curl_destination", "unitsquare_Regge_source", # This test fails in complex mode pytest.param("spheresphere", marks=pytest.mark.skipcomplex), @@ -180,14 +180,14 @@ def parameters(request): V_src = FunctionSpace(m_src, "N1curl", 2) # Not point evaluation nodes V_dest = VectorFunctionSpace(m_dest, "CG", 4) V_dest_2 = VectorFunctionSpace(m_dest, "DQ", 2) - elif request.param == "unitsquare_SminusDiv_destination": - m_src, m_dest, coords = unitsquaresetup() + elif request.param == "unitsquare_RT_N1curl_destination": + m_src, m_dest, coords = unitsquaresetup(dest_quad=False) expr_src = 2 * SpatialCoordinate(m_src) expr_dest = 2 * SpatialCoordinate(m_dest) expected = 2 * coords V_src = VectorFunctionSpace(m_src, "CG", 2) V_dest = FunctionSpace(m_dest, "RT", 2) # Not point evaluation nodes - V_dest_2 = FunctionSpace(m_dest, "N1Curl", 2) # Not point evaluation nodes + V_dest_2 = FunctionSpace(m_dest, "N1curl", 2) # Not point evaluation nodes elif request.param == "unitsquare_Regge_source": m_src, m_dest, coords = unitsquaresetup() expr_src = outer(SpatialCoordinate(m_src), SpatialCoordinate(m_src)) From 5052e41a351ae30bd601be0c93cfad52e76c3f6e Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 27 Jan 2026 13:01:49 +0000 Subject: [PATCH 14/34] add zero-form test --- .../regression/test_cross_mesh_non_lagrange.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 4796cb2452..85485fd3a0 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -96,7 +96,7 @@ def test_cross_mesh(V, rank): assert np.allclose(f_interpolated_direct.dat.data_ro, f_direct.dat.data_ro) -@pytest.mark.parametrize("rank", [1, 2]) +@pytest.mark.parametrize("rank", [0, 1, 2]) def test_cross_mesh_adjoint(V, rank): # Can already do Lagrange -> RT adjoint # V^* -> Q^* -> V_target^* @@ -117,7 +117,7 @@ def test_cross_mesh_adjoint(V, rank): expr = x2 target_expr = x1 - oneform_V = inner(expr, TestFunction(V)) * dx + oneform_V = inner(expr, TestFunction(V)) * dx # V^* cofunc_Vtarget_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) Q = make_quadrature_space(V) @@ -148,3 +148,10 @@ def test_cross_mesh_adjoint(V, rank): cofunc_Vtarget = assemble(interpolate(TestFunction(V_target), oneform_V)) # V^* -> V_target^* assert np.allclose(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) + elif rank == 0: + res = assemble(interpolate(target_expr, oneform_V)) + actual = assemble(inner(expr, expr) * dx) + assert np.isclose(res, actual) + +if __name__ == "__main__": + pytest.main([__file__ + "::test_cross_mesh_adjoint[RT_2-1]"]) From abfbb0c1d64a57ef78a24510d9d877db402e0946 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 27 Jan 2026 13:57:22 +0000 Subject: [PATCH 15/34] fix tests --- .../test_cross_mesh_non_lagrange.py | 30 +++++++++++++++---- 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 85485fd3a0..5153e36dee 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -5,6 +5,13 @@ from functools import partial +def mat_equals(a, b) -> bool: + """Check that two Matrices are equal.""" + a = a.petscmat.copy() + a.axpy(-1.0, b.petscmat) + return a.norm(norm_type=PETSc.NormType.NORM_FROBENIUS) < 1e-14 + + def fs_shape(V): shape = V.ufl_function_space().value_shape if len(shape) == 1: @@ -31,9 +38,9 @@ def make_quadrature_space(V): return fs_shape(V)(V.mesh(), element) -@pytest.fixture(params=[("RT", 2), ("RT", 3), ("RT", 4), ("BDM", 1), ("BDM", 2), ("BDM", 3), - ("BDFM", 2), ("HHJ", 2), ("N1curl", 2), ("N1curl", 3), ("N1curl", 4), - ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 2), ("GLS", 3), +@pytest.fixture(params=[("RT", 1), ("RT", 2), ("RT", 3), ("RT", 4), ("BDM", 1), ("BDM", 2), ("BDM", 3), + ("BDFM", 2), ("HHJ", 2), ("N1curl", 1), ("N1curl", 2), ("N1curl", 3), ("N1curl", 4), + ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 1), ("GLS", 2), ("GLS", 3), ("GLS", 4), ("GLS2", 1), ("GLS2", 2), ("GLS2", 3)], ids=lambda x: f"{x[0]}_{x[1]}") def V(request): @@ -42,6 +49,7 @@ def V(request): return FunctionSpace(mesh, element, degree) +@pytest.mark.parallel([1, 3]) @pytest.mark.parametrize("rank", [1, 2]) def test_cross_mesh(V, rank): mesh1 = UnitSquareMesh(5, 5) @@ -76,6 +84,7 @@ def test_cross_mesh(V, rank): # Direct assembly I_direct = assemble(interpolate(TrialFunction(V_source), V)) # V_source assert I_direct.arguments() == (TestFunction(V.dual()), TrialFunction(V_source)) + assert mat_equals(I_manual, I_direct) f_interpolated_manual = assemble(action(I_manual, f_source)) assert np.allclose(f_interpolated_manual.dat.data_ro, f_direct.dat.data_ro) @@ -96,10 +105,19 @@ def test_cross_mesh(V, rank): assert np.allclose(f_interpolated_direct.dat.data_ro, f_direct.dat.data_ro) +@pytest.mark.parallel([1, 3]) @pytest.mark.parametrize("rank", [0, 1, 2]) def test_cross_mesh_adjoint(V, rank): # Can already do Lagrange -> RT adjoint # V^* -> Q^* -> V_target^* + name = V.ufl_element()._short_name + deg = V.ufl_element().degree() + if name in ["N1curl", "GLS", "RT"] and deg == 1: + if name == "RT" and rank == 0: + pass + else: + pytest.skip(f"Not exact for degree {deg} {name} elements") + mesh1 = UnitSquareMesh(2, 2) x1 = SpatialCoordinate(mesh1) V_target = fs_shape(V)(mesh1, "CG", 1) @@ -131,6 +149,7 @@ def test_cross_mesh_adjoint(V, rank): # Direct assembly I_direct = assemble(interpolate(TestFunction(V_target), V)) # V^* x V_target -> R assert I_direct.arguments() == (TestFunction(V_target), TrialFunction(V.dual())) + assert mat_equals(I_manual, I_direct) cofunc_Vtarget_manual = assemble(action(I_manual, oneform_V)) assert np.allclose(cofunc_Vtarget_manual.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) @@ -149,9 +168,8 @@ def test_cross_mesh_adjoint(V, rank): cofunc_Vtarget = assemble(interpolate(TestFunction(V_target), oneform_V)) # V^* -> V_target^* assert np.allclose(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) elif rank == 0: + if name == "GLS2" and deg == 1: + pytest.skip(f"Not exact for degree {deg} {name} elements") res = assemble(interpolate(target_expr, oneform_V)) actual = assemble(inner(expr, expr) * dx) assert np.isclose(res, actual) - -if __name__ == "__main__": - pytest.main([__file__ + "::test_cross_mesh_adjoint[RT_2-1]"]) From 1ab45b993fca96ba97315f4e8748b7270d0b31ba Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 27 Jan 2026 13:57:26 +0000 Subject: [PATCH 16/34] docs --- docs/source/interpolation.rst | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/docs/source/interpolation.rst b/docs/source/interpolation.rst index 116b34cbef..985b0a258f 100644 --- a/docs/source/interpolation.rst +++ b/docs/source/interpolation.rst @@ -267,8 +267,16 @@ function space has finite elements (as given in the list of * **Lagrange/CG** (also known as Continuous Galerkin or P elements), * **Q** (i.e. Lagrange/CG on lines, quadrilaterals and hexahedra), -* **Discontinuous Lagrange/DG** (also known as Discontinuous Galerkin or DP elements) and -* **DQ** (i.e. Discontinuous Lagrange/DG on lines, quadrilaterals and hexahedra). +* **Discontinuous Lagrange/DG** (also known as Discontinuous Galerkin or DP elements), +* **DQ** (i.e. Discontinuous Lagrange/DG on lines, quadrilaterals and hexahedra), +* **Raviart-Thomas/RT**, +* **Brezzi-Douglas-Marini/BDM**, +* **Brezzi-Douglas-Fortin-Marini/BDFM**, +* **Hellan-Herrmann-Johnson/HHJ**, +* **Nedelec 1st kind/N1curl**, +* **Nedelec 2nd kind/N2curl**, +* **Gopalakrishnan-Lederer-Schoberl 1st kind/GLS**, +* **Gopalakrishnan-Lederer-Schoberl 2nd kind/GLS2**. Vector, tensor, and mixed function spaces can also be interpolated into from other meshes as long as they are constructed from these spaces. From 78b857e947bc7621e504d9e986c0e6ddc3c73ab6 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 27 Jan 2026 15:15:35 +0000 Subject: [PATCH 17/34] fixes --- firedrake/interpolation.py | 50 +++++++++---------- .../regression/test_interpolate_cross_mesh.py | 25 ---------- 2 files changed, 23 insertions(+), 52 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index c853ac326f..4821a7c320 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -530,7 +530,7 @@ def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpo If the target space does not have point-evaluation dofs. """ from firedrake.assemble import assemble - if target_space.ufl_element().mapping() != "identity": + if target_space.finat_element.mapping != "affine": raise ValueError(f"FunctionSpace {target_space} must have point-evaluation dofs.") # Immerse coordinates of target space point evaluation dofs in src_mesh @@ -571,18 +571,12 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) mat_type = mat_type or "aij" # Interpolate into intermediate quadrature space for non-identity mapped elements - if into_quadrature_space := self.target_space.ufl_element().mapping() != "identity": - Q_dest = self._get_quadrature_space(self.target_space) + if into_quadrature_space := self.target_space.finat_element.mapping != "affine": + target_space = self._get_quadrature_space(self.target_space) + f = Function(target_space.dual() if self.ufl_interpolate.is_adjoint else target_space) else: - Q_dest = None - target_space = Q_dest or self.target_space - - if into_quadrature_space and not self.ufl_interpolate.is_adjoint: - f = Function(target_space) - else: - # self.ufl_interpolate.function_space() is None in the 0-form case - f = Function(self.ufl_interpolate.function_space() or target_space) - # f = tensor or Function(target_space.dual() if self.ufl_interpolate.is_adjoint else target_space) + target_space = self.target_space + f = tensor or Function(self.ufl_interpolate.function_space() or target_space) point_eval, point_eval_input_ordering = self._get_symbolic_expressions(target_space) P0DG_vom_input_ordering = point_eval_input_ordering.argument_slots()[0].function_space().dual() @@ -600,12 +594,13 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) def callable() -> PETSc.Mat: res = assemble(interp_expr, mat_type=mat_type).petscmat if into_quadrature_space: + source_space = self.ufl_interpolate.function_space() if self.ufl_interpolate.is_adjoint: - I = AssembledMatrix((Argument(Q_dest, 0), Argument(target_space.dual(), 1)), None, res) - return assemble(action(I, interpolate(TestFunction(Q_dest), self.target_space))).petscmat + I = AssembledMatrix((Argument(source_space, 0), Argument(target_space.dual(), 1)), None, res) + return assemble(action(I, interpolate(TestFunction(target_space), self.target_space))).petscmat else: - I = AssembledMatrix((Argument(Q_dest.dual(), 0), Argument(target_space, 1)), None, res) - return assemble(action(interpolate(TrialFunction(Q_dest), self.target_space), I)).petscmat + I = AssembledMatrix((Argument(target_space.dual(), 0), Argument(source_space, 1)), None, res) + return assemble(action(interpolate(TrialFunction(target_space), self.target_space), I)).petscmat else: return res @@ -637,19 +632,20 @@ def callable() -> Cofunction: return f_target else: assert self.rank in {0, 1} - # We create the input-ordering Function before interpolating so we can - # set default missing values if required. - f_point_eval_input_ordering = Function(P0DG_vom_input_ordering) - if self.default_missing_val is not None: - f_point_eval_input_ordering.assign(self.default_missing_val) - elif self.allow_missing_dofs: - # If we allow missing points there may be points in the target - # mesh that are not in the source mesh. If we don't specify a - # default missing value we set these to NaN so we can identify - # them later. - f_point_eval_input_ordering.dat.data_wo[:] = numpy.nan def callable() -> Function | Number: + # We create the input-ordering Function before interpolating so we can + # set default missing values if required. + f_point_eval_input_ordering = Function(P0DG_vom_input_ordering) + if self.default_missing_val is not None: + f_point_eval_input_ordering.assign(self.default_missing_val) + elif self.allow_missing_dofs: + # If we allow missing points there may be points in the target + # mesh that are not in the source mesh. If we don't specify a + # default missing value we set these to NaN so we can identify + # them later. + f_point_eval_input_ordering.dat.data_wo[:] = numpy.nan + assemble(action(point_eval_input_ordering, point_eval), tensor=f_point_eval_input_ordering) # We assign these values to the output function if self.allow_missing_dofs and self.default_missing_val is None: diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index b721fb40d2..61f4992577 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -410,31 +410,6 @@ def test_interpolate_unitsquare_tfs_shape(shape, symmetry): assemble(interpolate(f_src, V_dest)) -def test_interpolate_cross_mesh_not_point_eval(): - m_src = UnitSquareMesh(2, 3) - m_dest = UnitSquareMesh(3, 5, quadrilateral=True) - coords = np.array( - [[0.56, 0.6], [0.1, 0.9], [0.9, 0.1], [0.9, 0.9], [0.726, 0.6584]] - ) # fairly arbitrary - # add the coordinates of the mesh vertices to test boundaries - vertices_src = allgather(m_src.comm, m_src.coordinates.dat.data_ro) - coords = np.concatenate((coords, vertices_src)) - vertices_dest = allgather(m_dest.comm, m_dest.coordinates.dat.data_ro) - coords = np.concatenate((coords, vertices_dest)) - dest_eval = PointEvaluator(m_dest, coords) - expr_src = 2 * SpatialCoordinate(m_src) - expr_dest = 2 * SpatialCoordinate(m_dest) - expected = 2 * coords - V_src = FunctionSpace(m_src, "RT", 2) - V_dest = FunctionSpace(m_dest, "RTCE", 2) - atol = 1e-8 # default - # This might not make much mathematical sense, but it should test if we get - # the not implemented error for non-point evaluation nodes! - interpolate_function( - m_src, m_dest, V_src, V_dest, dest_eval, expected, expr_src, expr_dest, atol - ) - - def interpolate_function( m_src, m_dest, V_src, V_dest, dest_eval, expected, expr_src, expr_dest, atol ): From 73ede41af748403199859d1d67beadf33eb7f84d Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 28 Jan 2026 13:51:16 +0000 Subject: [PATCH 18/34] fix test --- tests/firedrake/regression/test_interpolate_cross_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index 61f4992577..a04e3edbc0 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -39,7 +39,7 @@ def make_high_order(m_low_order, degree): @pytest.fixture( params=[ - "unitsquare_RT_N1curl_destination", + "unitsquare", "circlemanifold", "circlemanifold_to_high_order", "unitsquare_from_high_order", From 2ad7fcf139db6b9c2bb7dad50d12b70d03b0ab8a Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Thu, 29 Jan 2026 14:59:33 +0000 Subject: [PATCH 19/34] review suggestions --- firedrake/interpolation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 4821a7c320..378042801d 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -530,7 +530,7 @@ def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpo If the target space does not have point-evaluation dofs. """ from firedrake.assemble import assemble - if target_space.finat_element.mapping != "affine": + if not target_space.finat_element.has_pointwise_dual_basis: raise ValueError(f"FunctionSpace {target_space} must have point-evaluation dofs.") # Immerse coordinates of target space point evaluation dofs in src_mesh @@ -571,7 +571,7 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) mat_type = mat_type or "aij" # Interpolate into intermediate quadrature space for non-identity mapped elements - if into_quadrature_space := self.target_space.finat_element.mapping != "affine": + if into_quadrature_space := not self.target_space.finat_element.has_pointwise_dual_basis: target_space = self._get_quadrature_space(self.target_space) f = Function(target_space.dual() if self.ufl_interpolate.is_adjoint else target_space) else: From e22e2076b9822d94e0815d830dbcf091aafdf564 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Thu, 29 Jan 2026 18:46:56 +0000 Subject: [PATCH 20/34] automate docs --- docs/source/element_list.py | 4 +++- docs/source/variational-problems.rst | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/source/element_list.py b/docs/source/element_list.py index 4671bfcf1b..9c7bc50af7 100644 --- a/docs/source/element_list.py +++ b/docs/source/element_list.py @@ -32,4 +32,6 @@ def cells(cellnames): cellnames = cells(cellnames) shape = shape_names[value_rank] - csvwriter.writerow((family, short_name, shape, cellnames)) + interpolatable = hasattr(element, "dual_basis") + + csvwriter.writerow((family, short_name, shape, cellnames, interpolatable)) diff --git a/docs/source/variational-problems.rst b/docs/source/variational-problems.rst index e164f6c620..2f9ee9cddc 100644 --- a/docs/source/variational-problems.rst +++ b/docs/source/variational-problems.rst @@ -270,7 +270,7 @@ Supported finite elements Firedrake supports the use of the following finite elements. .. csv-table:: - :header: "Name", "Short name", "Value shape", "Valid cells" + :header: "Name", "Short name", "Value shape", "Valid cells", "Can interpolate into?" :widths: 20, 10, 10, 40 :file: element_list.csv From 7eb6e6934b9cd988a0b8c0906c534d02a66ad163 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Fri, 30 Jan 2026 16:40:35 +0000 Subject: [PATCH 21/34] review suggestions fixes fixes --- docs/source/element_list.py | 4 +- docs/source/interpolation.rst | 2 +- docs/source/variational-problems.rst | 2 +- firedrake/__init__.py | 2 +- firedrake/exceptions.py | 10 ++++- firedrake/interpolation.py | 44 ++++--------------- .../test_cross_mesh_non_lagrange.py | 21 +++++---- .../regression/test_interpolate_cross_mesh.py | 2 +- .../regression/test_interpolation_manual.py | 8 ++-- 9 files changed, 38 insertions(+), 57 deletions(-) diff --git a/docs/source/element_list.py b/docs/source/element_list.py index 9c7bc50af7..4671bfcf1b 100644 --- a/docs/source/element_list.py +++ b/docs/source/element_list.py @@ -32,6 +32,4 @@ def cells(cellnames): cellnames = cells(cellnames) shape = shape_names[value_rank] - interpolatable = hasattr(element, "dual_basis") - - csvwriter.writerow((family, short_name, shape, cellnames, interpolatable)) + csvwriter.writerow((family, short_name, shape, cellnames)) diff --git a/docs/source/interpolation.rst b/docs/source/interpolation.rst index 985b0a258f..94f9290466 100644 --- a/docs/source/interpolation.rst +++ b/docs/source/interpolation.rst @@ -325,7 +325,7 @@ Interpolating onto other meshes ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ If the target mesh extends outside the source mesh domain, then cross-mesh -interpolation will raise a :py:class:`~.DofNotDefinedError`. +interpolation will raise a :py:class:`~.DoFNotDefinedError`. .. literalinclude:: ../../tests/firedrake/regression/test_interpolation_manual.py :language: python3 diff --git a/docs/source/variational-problems.rst b/docs/source/variational-problems.rst index 2f9ee9cddc..e164f6c620 100644 --- a/docs/source/variational-problems.rst +++ b/docs/source/variational-problems.rst @@ -270,7 +270,7 @@ Supported finite elements Firedrake supports the use of the following finite elements. .. csv-table:: - :header: "Name", "Short name", "Value shape", "Valid cells", "Can interpolate into?" + :header: "Name", "Short name", "Value shape", "Valid cells" :widths: 20, 10, 10, 40 :file: element_list.csv diff --git a/firedrake/__init__.py b/firedrake/__init__.py index 1238563ba4..300c32880e 100644 --- a/firedrake/__init__.py +++ b/firedrake/__init__.py @@ -69,7 +69,7 @@ def init_petsc(): from firedrake.deflation import DeflatedSNES, Deflation # noqa: F401 from firedrake.exceptions import ( # noqa: F401 FiredrakeException, ConvergenceError, MismatchingDomainError, - VertexOnlyMeshMissingPointsError, DofNotDefinedError + VertexOnlyMeshMissingPointsError, DoFNotDefinedError, DoFTypeError, ) from firedrake.function import ( # noqa: F401 Function, PointNotInDomainError, diff --git a/firedrake/exceptions.py b/firedrake/exceptions.py index 8816375138..fb93686d5c 100644 --- a/firedrake/exceptions.py +++ b/firedrake/exceptions.py @@ -9,14 +9,20 @@ class ConvergenceError(FiredrakeException): """Error raised when a solver fails to converge.""" -class DofNotDefinedError(FiredrakeException): +class DoFNotDefinedError(FiredrakeException): r"""Raised when attempting to interpolate across function spaces where the - target function space contains degrees of freedom (i.e. nodes) which cannot + target function space contains degrees of freedom (i.e. DoFs) which cannot be defined in the source function space. This typically occurs when the target mesh covers a larger domain than the source mesh. """ +class DoFTypeError(FiredrakeException): + """Raised when an operation is attempted on a degree of freedom (DoF) + type which is not supported. + """ + + class VertexOnlyMeshMissingPointsError(FiredrakeException): """Exception raised when 1 or more points are not found by a :func:`~.VertexOnlyMesh` in its parent mesh. diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 378042801d..1fd5be987b 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -41,7 +41,8 @@ from firedrake.function import Function from firedrake.cofunction import Cofunction from firedrake.exceptions import ( - DofNotDefinedError, VertexOnlyMeshMissingPointsError, NonUniqueMeshSequenceError + DoFNotDefinedError, VertexOnlyMeshMissingPointsError, NonUniqueMeshSequenceError, + DoFTypeError, ) from mpi4py import MPI @@ -448,13 +449,7 @@ def _get_element(self, V: WithGeometry) -> FiniteElementBase: if isinstance(dest_element, MixedElement): if isinstance(dest_element, VectorElement | TensorElement): # In this case all sub elements are equal - base_element = dest_element.sub_elements[0] - if base_element.reference_value_shape != (): - raise NotImplementedError( - "Can't yet cross-mesh interpolate onto function spaces made from VectorElements " - "or TensorElements made from sub elements with value shape other than ()." - ) - return base_element + return dest_element.sub_elements[0] else: raise NotImplementedError("Interpolation with MixedFunctionSpace requires MixedInterpolator.") else: @@ -484,29 +479,6 @@ def _fs_type(self, V: WithGeometry) -> Callable[..., WithGeometry]: symmetry = V.ufl_element().symmetry() return partial(TensorFunctionSpace, shape=shape, symmetry=symmetry) - def _get_quadrature_space(self, V: WithGeometry) -> WithGeometry: - """Return a FunctionSpace whose element is a quadrature element with - points at the quadrature points of V's element. - - Used as an intermediate space for interpolating into a space which doesn't - have point evaluation dofs. - - Parameters - ---------- - V : WithGeometry - A :class:`.WithGeometry` function space to build the quadrature space for. - - Returns - ------- - WithGeometry - A :class:`.WithGeometry` function space with quadrature element. - """ - V_element = V.finat_element - _, points = V_element.dual_basis - weights = numpy.full(len(points.points), numpy.nan) # These can be any number since we never integrate - quad_scheme = QuadratureRule(points, weights, ref_el=V_element.cell) - return self._fs_type(V)(V.mesh(), "Quadrature", degree=V_element.degree, quad_scheme=quad_scheme) - def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpolate, Interpolate]: """Return symbolic ``Interpolate`` expressions for point evaluation of the `target_space`s dofs in the source mesh, and the corresponding input-ordering interpolation. @@ -524,14 +496,14 @@ def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpo Raises ------ - DofNotDefinedError + DoFNotDefinedError If any of the target spaces dofs cannot be defined in the source mesh. - ValueError + DoFTypeError If the target space does not have point-evaluation dofs. """ from firedrake.assemble import assemble if not target_space.finat_element.has_pointwise_dual_basis: - raise ValueError(f"FunctionSpace {target_space} must have point-evaluation dofs.") + raise DoFTypeError(f"FunctionSpace {target_space} must have point-evaluation dofs.") # Immerse coordinates of target space point evaluation dofs in src_mesh target_mesh = target_space.mesh().unique() @@ -546,7 +518,7 @@ def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpo missing_points_behaviour=self.missing_points_behaviour, ) except VertexOnlyMeshMissingPointsError: - raise DofNotDefinedError(f"The given target function space on domain {self.target_mesh} " + raise DoFNotDefinedError(f"The given target function space on domain {self.target_mesh} " "contains degrees of freedom which cannot cannot be defined in the " f"source function space on domain {self.source_mesh}. " "This may be because the target mesh covers a larger domain than the " @@ -572,7 +544,7 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) # Interpolate into intermediate quadrature space for non-identity mapped elements if into_quadrature_space := not self.target_space.finat_element.has_pointwise_dual_basis: - target_space = self._get_quadrature_space(self.target_space) + target_space = self.target_space.quadrature_space() f = Function(target_space.dual() if self.ufl_interpolate.is_adjoint else target_space) else: target_space = self.target_space diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 5153e36dee..e858b69a6d 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -14,13 +14,14 @@ def mat_equals(a, b) -> bool: def fs_shape(V): shape = V.ufl_function_space().value_shape - if len(shape) == 1: - fs_type = partial(VectorFunctionSpace, dim=shape[0]) + if len(shape) == 0: + return FunctionSpace + elif len(shape) == 1: + return partial(VectorFunctionSpace, dim=shape[0]) elif len(shape) == 2: - fs_type = partial(TensorFunctionSpace, shape=shape) + return partial(TensorFunctionSpace, shape=shape) else: raise ValueError("Invalid function space shape") - return fs_type def make_quadrature_space(V): @@ -38,7 +39,7 @@ def make_quadrature_space(V): return fs_shape(V)(V.mesh(), element) -@pytest.fixture(params=[("RT", 1), ("RT", 2), ("RT", 3), ("RT", 4), ("BDM", 1), ("BDM", 2), ("BDM", 3), +@pytest.fixture(params=[("ARG", 5), ("RT", 1), ("RT", 2), ("RT", 3), ("RT", 4), ("BDM", 1), ("BDM", 2), ("BDM", 3), ("BDFM", 2), ("HHJ", 2), ("N1curl", 1), ("N1curl", 2), ("N1curl", 3), ("N1curl", 4), ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 1), ("GLS", 2), ("GLS", 3), ("GLS", 4), ("GLS2", 1), ("GLS2", 2), ("GLS2", 3)], @@ -49,7 +50,7 @@ def V(request): return FunctionSpace(mesh, element, degree) -@pytest.mark.parallel([1, 3]) +# @pytest.mark.parallel([1, 3]) @pytest.mark.parametrize("rank", [1, 2]) def test_cross_mesh(V, rank): mesh1 = UnitSquareMesh(5, 5) @@ -58,7 +59,11 @@ def test_cross_mesh(V, rank): x1, y1 = SpatialCoordinate(mesh2) shape = V.ufl_function_space().value_shape - if len(shape) == 1: + if len(shape) == 0: + fs_type = FunctionSpace + expr1 = x * x + y * y + expr2 = x1 * x1 + y1 * y1 + elif len(shape) == 1: fs_type = partial(VectorFunctionSpace, dim=shape[0]) expr1 = as_vector([x, y]) expr2 = as_vector([x1, y1]) @@ -105,7 +110,7 @@ def test_cross_mesh(V, rank): assert np.allclose(f_interpolated_direct.dat.data_ro, f_direct.dat.data_ro) -@pytest.mark.parallel([1, 3]) +# @pytest.mark.parallel([1, 3]) @pytest.mark.parametrize("rank", [0, 1, 2]) def test_cross_mesh_adjoint(V, rank): # Can already do Lagrange -> RT adjoint diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index a04e3edbc0..94cb4a81b3 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -528,7 +528,7 @@ def test_missing_dofs(): expr = x * y V_src = FunctionSpace(m_src, "CG", 2) V_dest = FunctionSpace(m_dest, "CG", 3) - with pytest.raises(DofNotDefinedError): + with pytest.raises(DoFNotDefinedError): assemble(interpolate(TrialFunction(V_src), V_dest)) f_src = Function(V_src).interpolate(expr) f_dest = assemble(interpolate(f_src, V_dest, allow_missing_dofs=True)) diff --git a/tests/firedrake/regression/test_interpolation_manual.py b/tests/firedrake/regression/test_interpolation_manual.py index d6dfcf09d6..10704ef359 100644 --- a/tests/firedrake/regression/test_interpolation_manual.py +++ b/tests/firedrake/regression/test_interpolation_manual.py @@ -158,13 +158,13 @@ def correct_indent(): src_mesh, dest_mesh, f_src, V_dest = correct_indent() - with pytest.raises(DofNotDefinedError): + with pytest.raises(DoFNotDefinedError): # [test_cross_mesh 3] - # ... but get a DofNotDefinedError if we try - f_dest = assemble(interpolate(f_src, V_dest)) # raises DofNotDefinedError + # ... but get a DoFNotDefinedError if we try + f_dest = assemble(interpolate(f_src, V_dest)) # raises DoFNotDefinedError # [test_cross_mesh 4] - with pytest.raises(DofNotDefinedError): + with pytest.raises(DoFNotDefinedError): # as will the interpolate method of a Function f_dest = Function(V_dest).interpolate(f_src) From f7ecb15cf62805b6086f347f296cd04d29977202 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 2 Feb 2026 13:36:58 +0000 Subject: [PATCH 22/34] fix --- tests/firedrake/regression/test_cross_mesh_non_lagrange.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index e858b69a6d..d33133631c 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -39,7 +39,7 @@ def make_quadrature_space(V): return fs_shape(V)(V.mesh(), element) -@pytest.fixture(params=[("ARG", 5), ("RT", 1), ("RT", 2), ("RT", 3), ("RT", 4), ("BDM", 1), ("BDM", 2), ("BDM", 3), +@pytest.fixture(params=[("RT", 1), ("RT", 2), ("RT", 3), ("RT", 4), ("BDM", 1), ("BDM", 2), ("BDM", 3), ("BDFM", 2), ("HHJ", 2), ("N1curl", 1), ("N1curl", 2), ("N1curl", 3), ("N1curl", 4), ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 1), ("GLS", 2), ("GLS", 3), ("GLS", 4), ("GLS2", 1), ("GLS2", 2), ("GLS2", 3)], From 01d65facb50ab700d50c8f1a7f3bf7f7d551cff9 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 2 Feb 2026 15:51:28 +0000 Subject: [PATCH 23/34] undo API changes fixe --- docs/source/interpolation.rst | 14 +++----------- firedrake/__init__.py | 2 +- firedrake/exceptions.py | 6 +++--- firedrake/interpolation.py | 8 ++++---- .../regression/test_interpolate_cross_mesh.py | 2 +- .../regression/test_interpolation_manual.py | 13 +++++++------ 6 files changed, 19 insertions(+), 26 deletions(-) diff --git a/docs/source/interpolation.rst b/docs/source/interpolation.rst index 94f9290466..116b34cbef 100644 --- a/docs/source/interpolation.rst +++ b/docs/source/interpolation.rst @@ -267,16 +267,8 @@ function space has finite elements (as given in the list of * **Lagrange/CG** (also known as Continuous Galerkin or P elements), * **Q** (i.e. Lagrange/CG on lines, quadrilaterals and hexahedra), -* **Discontinuous Lagrange/DG** (also known as Discontinuous Galerkin or DP elements), -* **DQ** (i.e. Discontinuous Lagrange/DG on lines, quadrilaterals and hexahedra), -* **Raviart-Thomas/RT**, -* **Brezzi-Douglas-Marini/BDM**, -* **Brezzi-Douglas-Fortin-Marini/BDFM**, -* **Hellan-Herrmann-Johnson/HHJ**, -* **Nedelec 1st kind/N1curl**, -* **Nedelec 2nd kind/N2curl**, -* **Gopalakrishnan-Lederer-Schoberl 1st kind/GLS**, -* **Gopalakrishnan-Lederer-Schoberl 2nd kind/GLS2**. +* **Discontinuous Lagrange/DG** (also known as Discontinuous Galerkin or DP elements) and +* **DQ** (i.e. Discontinuous Lagrange/DG on lines, quadrilaterals and hexahedra). Vector, tensor, and mixed function spaces can also be interpolated into from other meshes as long as they are constructed from these spaces. @@ -325,7 +317,7 @@ Interpolating onto other meshes ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ If the target mesh extends outside the source mesh domain, then cross-mesh -interpolation will raise a :py:class:`~.DoFNotDefinedError`. +interpolation will raise a :py:class:`~.DofNotDefinedError`. .. literalinclude:: ../../tests/firedrake/regression/test_interpolation_manual.py :language: python3 diff --git a/firedrake/__init__.py b/firedrake/__init__.py index 300c32880e..de15331684 100644 --- a/firedrake/__init__.py +++ b/firedrake/__init__.py @@ -69,7 +69,7 @@ def init_petsc(): from firedrake.deflation import DeflatedSNES, Deflation # noqa: F401 from firedrake.exceptions import ( # noqa: F401 FiredrakeException, ConvergenceError, MismatchingDomainError, - VertexOnlyMeshMissingPointsError, DoFNotDefinedError, DoFTypeError, + VertexOnlyMeshMissingPointsError, DofNotDefinedError, DofTypeError, ) from firedrake.function import ( # noqa: F401 Function, PointNotInDomainError, diff --git a/firedrake/exceptions.py b/firedrake/exceptions.py index fb93686d5c..1a74c6a35e 100644 --- a/firedrake/exceptions.py +++ b/firedrake/exceptions.py @@ -9,15 +9,15 @@ class ConvergenceError(FiredrakeException): """Error raised when a solver fails to converge.""" -class DoFNotDefinedError(FiredrakeException): +class DofNotDefinedError(FiredrakeException): r"""Raised when attempting to interpolate across function spaces where the - target function space contains degrees of freedom (i.e. DoFs) which cannot + target function space contains degrees of freedom (i.e. nodes) which cannot be defined in the source function space. This typically occurs when the target mesh covers a larger domain than the source mesh. """ -class DoFTypeError(FiredrakeException): +class DofTypeError(FiredrakeException): """Raised when an operation is attempted on a degree of freedom (DoF) type which is not supported. """ diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 1fd5be987b..28c75dfdbd 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -41,8 +41,8 @@ from firedrake.function import Function from firedrake.cofunction import Cofunction from firedrake.exceptions import ( - DoFNotDefinedError, VertexOnlyMeshMissingPointsError, NonUniqueMeshSequenceError, - DoFTypeError, + DofNotDefinedError, VertexOnlyMeshMissingPointsError, NonUniqueMeshSequenceError, + DofTypeError, ) from mpi4py import MPI @@ -503,7 +503,7 @@ def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpo """ from firedrake.assemble import assemble if not target_space.finat_element.has_pointwise_dual_basis: - raise DoFTypeError(f"FunctionSpace {target_space} must have point-evaluation dofs.") + raise DofTypeError(f"FunctionSpace {target_space} must have point-evaluation dofs.") # Immerse coordinates of target space point evaluation dofs in src_mesh target_mesh = target_space.mesh().unique() @@ -518,7 +518,7 @@ def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpo missing_points_behaviour=self.missing_points_behaviour, ) except VertexOnlyMeshMissingPointsError: - raise DoFNotDefinedError(f"The given target function space on domain {self.target_mesh} " + raise DofNotDefinedError(f"The given target function space on domain {self.target_mesh} " "contains degrees of freedom which cannot cannot be defined in the " f"source function space on domain {self.source_mesh}. " "This may be because the target mesh covers a larger domain than the " diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index 94cb4a81b3..a04e3edbc0 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -528,7 +528,7 @@ def test_missing_dofs(): expr = x * y V_src = FunctionSpace(m_src, "CG", 2) V_dest = FunctionSpace(m_dest, "CG", 3) - with pytest.raises(DoFNotDefinedError): + with pytest.raises(DofNotDefinedError): assemble(interpolate(TrialFunction(V_src), V_dest)) f_src = Function(V_src).interpolate(expr) f_dest = assemble(interpolate(f_src, V_dest, allow_missing_dofs=True)) diff --git a/tests/firedrake/regression/test_interpolation_manual.py b/tests/firedrake/regression/test_interpolation_manual.py index 10704ef359..cedaee54be 100644 --- a/tests/firedrake/regression/test_interpolation_manual.py +++ b/tests/firedrake/regression/test_interpolation_manual.py @@ -102,6 +102,7 @@ def mydata(points): def test_line_integral(): # [test_line_integral 1] + from firedrake.mesh import Mesh, plex_from_cell_list # Start with a simple field exactly represented in the function space over # the unit square domain. m = UnitSquareMesh(2, 2) @@ -113,8 +114,8 @@ def test_line_integral(): # Note that it only has 1 cell cells = np.asarray([[0, 1]]) vertex_coords = np.asarray([[0.0, 0.0], [1.0, 1.0]]) - plex = mesh.plex_from_cell_list(1, cells, vertex_coords, comm=m.comm) - line = mesh.Mesh(plex, dim=2) + plex = plex_from_cell_list(1, cells, vertex_coords, comm=m.comm) + line = Mesh(plex, dim=2) # [test_line_integral 2] x, y = SpatialCoordinate(line) V_line = FunctionSpace(line, "CG", 2) @@ -158,13 +159,13 @@ def correct_indent(): src_mesh, dest_mesh, f_src, V_dest = correct_indent() - with pytest.raises(DoFNotDefinedError): + with pytest.raises(DofNotDefinedError): # [test_cross_mesh 3] - # ... but get a DoFNotDefinedError if we try - f_dest = assemble(interpolate(f_src, V_dest)) # raises DoFNotDefinedError + # ... but get a DofNotDefinedError if we try + f_dest = assemble(interpolate(f_src, V_dest)) # raises DofNotDefinedError # [test_cross_mesh 4] - with pytest.raises(DoFNotDefinedError): + with pytest.raises(DofNotDefinedError): # as will the interpolate method of a Function f_dest = Function(V_dest).interpolate(f_src) From 4321ad691c2f3bbfe97d77e66751ebe7e8eb303a Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Thu, 5 Feb 2026 12:30:50 +0000 Subject: [PATCH 24/34] update tests --- .../test_cross_mesh_non_lagrange.py | 31 +++++-------------- 1 file changed, 8 insertions(+), 23 deletions(-) diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index d33133631c..e0daba122e 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -24,25 +24,10 @@ def fs_shape(V): raise ValueError("Invalid function space shape") -def make_quadrature_space(V): - """Builds a quadrature space on the target mesh. - - This space has point evaluation dofs at the quadrature PointSet - of the target space's element - - """ - fe = V.finat_element - _, ps = fe.dual_basis - wts = np.full(len(ps.points), np.nan) # These can be any number since we never integrate - scheme = QuadratureRule(ps, wts, ref_el=fe.cell) - element = FiniteElement("Quadrature", degree=fe.degree, quad_scheme=scheme) - return fs_shape(V)(V.mesh(), element) - - -@pytest.fixture(params=[("RT", 1), ("RT", 2), ("RT", 3), ("RT", 4), ("BDM", 1), ("BDM", 2), ("BDM", 3), - ("BDFM", 2), ("HHJ", 2), ("N1curl", 1), ("N1curl", 2), ("N1curl", 3), ("N1curl", 4), - ("N2curl", 1), ("N2curl", 2), ("N2curl", 3), ("GLS", 1), ("GLS", 2), ("GLS", 3), - ("GLS", 4), ("GLS2", 1), ("GLS2", 2), ("GLS2", 3)], +@pytest.fixture(params=[("RT", 1), ("RT", 2), ("BDM", 1), ("BDM", 2), ("BDFM", 2), + ("HHJ", 0), ("HHJ", 2), ("N1curl", 1), ("N1curl", 2), + ("N2curl", 1), ("N2curl", 2), ("GLS", 1), ("GLS", 2), + ("GLS2", 2), ("Regge", 0), ("Regge", 2)], ids=lambda x: f"{x[0]}_{x[1]}") def V(request): element, degree = request.param @@ -78,7 +63,7 @@ def test_cross_mesh(V, rank): f_source = Function(V_source).interpolate(expr1) f_direct = Function(V).interpolate(expr2) - Q = make_quadrature_space(V) + Q = V.quadrature_space() if rank == 2: # Assemble the operator @@ -122,6 +107,8 @@ def test_cross_mesh_adjoint(V, rank): pass else: pytest.skip(f"Not exact for degree {deg} {name} elements") + elif name in ["Regge", "HHJ"] and deg == 0: + pytest.skip(f"Not exact for degree {deg} {name} elements") mesh1 = UnitSquareMesh(2, 2) x1 = SpatialCoordinate(mesh1) @@ -143,7 +130,7 @@ def test_cross_mesh_adjoint(V, rank): oneform_V = inner(expr, TestFunction(V)) * dx # V^* cofunc_Vtarget_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) - Q = make_quadrature_space(V) + Q = V.quadrature_space() if rank == 2: # Assemble the operator @@ -173,8 +160,6 @@ def test_cross_mesh_adjoint(V, rank): cofunc_Vtarget = assemble(interpolate(TestFunction(V_target), oneform_V)) # V^* -> V_target^* assert np.allclose(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) elif rank == 0: - if name == "GLS2" and deg == 1: - pytest.skip(f"Not exact for degree {deg} {name} elements") res = assemble(interpolate(target_expr, oneform_V)) actual = assemble(inner(expr, expr) * dx) assert np.isclose(res, actual) From 5b4c9f4e71465d90959ce8d3024ebb8543273455 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Thu, 5 Feb 2026 13:11:53 +0000 Subject: [PATCH 25/34] cache symbolic expressions for cross-mesh interpolation --- firedrake/interpolation.py | 88 +++++++++++++++++++------------------- 1 file changed, 45 insertions(+), 43 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 28c75dfdbd..10758a9477 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -430,9 +430,22 @@ def __init__(self, expr: Interpolate): if self.source_mesh.unique().geometric_dimension != self.target_mesh.unique().geometric_dimension: raise ValueError("Geometric dimensions of source and destination meshes must match.") + + # Interpolate into intermediate quadrature space for non-identity mapped elements + if into_quadrature_space := not self.target_space.finat_element.has_pointwise_dual_basis: + self.original_target_space = self.target_space + r"""The original target space for interpolation, as specified by the user. + This is only used if ``self.into_quadrature_space`` is ``True``.""" + self.target_space = self.target_space.quadrature_space() + r"""The target space for the cross-mesh interpolation. Must have point-evaluation dofs. + If ``self.original_target_space`` does not have point-evaluation dofs, then this is + an intermediate quadrature space.""" - def _get_element(self, V: WithGeometry) -> FiniteElementBase: - """Return the element of the function space V. If V is tensor/vector valued, + self.into_quadrature_space = into_quadrature_space + + @cached_property + def _target_space_element(self) -> FiniteElementBase: + """Return the element of `self.target_space`. If V is tensor/vector valued, return the base scalar element. Parameters @@ -445,7 +458,7 @@ def _get_element(self, V: WithGeometry) -> FiniteElementBase: FiniteElementBase The base element of V. """ - dest_element = V.ufl_element() + dest_element = self.target_space.ufl_element() if isinstance(dest_element, MixedElement): if isinstance(dest_element, VectorElement | TensorElement): # In this case all sub elements are equal @@ -456,38 +469,30 @@ def _get_element(self, V: WithGeometry) -> FiniteElementBase: # scalar fiat/finat element return dest_element - def _fs_type(self, V: WithGeometry) -> Callable[..., WithGeometry]: - """Returns a callable which returns a function space matching the type of V. - - Parameters - ---------- - V - A :class:`.WithGeometry` function space. + @cached_property + def _target_space_type(self) -> Callable[..., WithGeometry]: + """Returns a callable which returns a function space matching the type of `self.target_space`. Returns ------- Callable - A callable which returns a :class:`.WithGeometry` matching the type of V. + A callable which returns a :class:`.WithGeometry` matching the type of `self.target_space`. """ # Get the correct type of function space - shape = V.value_shape + shape = self.target_space.value_shape if len(shape) == 0: return FunctionSpace elif len(shape) == 1: return partial(VectorFunctionSpace, dim=shape[0]) else: - symmetry = V.ufl_element().symmetry() + symmetry = self.target_space.ufl_element().symmetry() return partial(TensorFunctionSpace, shape=shape, symmetry=symmetry) - def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpolate, Interpolate]: - """Return symbolic ``Interpolate`` expressions for point evaluation of the `target_space`s + @cached_property + def _symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: + """Return symbolic ``Interpolate`` expressions for point evaluation of `self.target_space`s dofs in the source mesh, and the corresponding input-ordering interpolation. - Parameters - ---------- - target_space - The :class:`.WithGeometry` function space which we are interpolating into. - Returns ------- tuple[Interpolate, Interpolate] @@ -502,12 +507,12 @@ def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpo If the target space does not have point-evaluation dofs. """ from firedrake.assemble import assemble - if not target_space.finat_element.has_pointwise_dual_basis: - raise DofTypeError(f"FunctionSpace {target_space} must have point-evaluation dofs.") + if not self.target_space.finat_element.has_pointwise_dual_basis: + raise DofTypeError(f"FunctionSpace {self.target_space} must have point-evaluation dofs.") # Immerse coordinates of target space point evaluation dofs in src_mesh - target_mesh = target_space.mesh().unique() - target_space_vec = VectorFunctionSpace(target_mesh, self._get_element(target_space)) + target_mesh = self.target_space.mesh().unique() + target_space_vec = VectorFunctionSpace(target_mesh, self._target_space_element) f_dest_node_coords = assemble(interpolate(target_mesh.coordinates, target_space_vec)) dest_node_coords = f_dest_node_coords.dat.data_ro.reshape(-1, target_mesh.geometric_dimension) try: @@ -525,12 +530,11 @@ def _get_symbolic_expressions(self, target_space: WithGeometry) -> tuple[Interpo "source mesh. To disable this error, set allow_missing_dofs=True.") # Get expression for point evaluation at the dest_node_coords - fs_type = self._fs_type(target_space) - P0DG_vom = fs_type(vom, "DG", 0) + P0DG_vom = self._target_space_type(vom, "DG", 0) point_eval = interpolate(self.operand, P0DG_vom) # Interpolate into the input-ordering VOM - P0DG_vom_input_ordering = fs_type(vom.input_ordering, "DG", 0) + P0DG_vom_input_ordering = self._target_space_type(vom.input_ordering, "DG", 0) arg = Argument(P0DG_vom, 0 if self.ufl_interpolate.is_adjoint else 1) point_eval_input_ordering = interpolate(arg, P0DG_vom_input_ordering) @@ -543,14 +547,12 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) mat_type = mat_type or "aij" # Interpolate into intermediate quadrature space for non-identity mapped elements - if into_quadrature_space := not self.target_space.finat_element.has_pointwise_dual_basis: - target_space = self.target_space.quadrature_space() - f = Function(target_space.dual() if self.ufl_interpolate.is_adjoint else target_space) + if self.into_quadrature_space: + f = Function(self.target_space.dual() if self.ufl_interpolate.is_adjoint else self.target_space) else: - target_space = self.target_space - f = tensor or Function(self.ufl_interpolate.function_space() or target_space) + f = tensor or Function(self.ufl_interpolate.function_space() or self.target_space) - point_eval, point_eval_input_ordering = self._get_symbolic_expressions(target_space) + point_eval, point_eval_input_ordering = self._symbolic_expressions P0DG_vom_input_ordering = point_eval_input_ordering.argument_slots()[0].function_space().dual() if self.rank == 2: @@ -565,14 +567,14 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) def callable() -> PETSc.Mat: res = assemble(interp_expr, mat_type=mat_type).petscmat - if into_quadrature_space: + if self.into_quadrature_space: source_space = self.ufl_interpolate.function_space() if self.ufl_interpolate.is_adjoint: - I = AssembledMatrix((Argument(source_space, 0), Argument(target_space.dual(), 1)), None, res) - return assemble(action(I, interpolate(TestFunction(target_space), self.target_space))).petscmat + I = AssembledMatrix((Argument(source_space, 0), Argument(self.target_space.dual(), 1)), None, res) + return assemble(action(I, interpolate(TestFunction(self.target_space), self.original_target_space))).petscmat else: - I = AssembledMatrix((Argument(target_space.dual(), 0), Argument(source_space, 1)), None, res) - return assemble(action(interpolate(TrialFunction(target_space), self.target_space), I)).petscmat + I = AssembledMatrix((Argument(self.target_space.dual(), 0), Argument(source_space, 1)), None, res) + return assemble(action(interpolate(TrialFunction(self.target_space), self.original_target_space), I)).petscmat else: return res @@ -580,8 +582,8 @@ def callable() -> PETSc.Mat: assert self.rank == 1 def callable() -> Cofunction: - if into_quadrature_space: - cofunc = assemble(interpolate(TestFunction(target_space), self.dual_arg)) + if self.into_quadrature_space: + cofunc = assemble(interpolate(TestFunction(self.target_space), self.dual_arg)) f_target = Cofunction(point_eval.function_space()) else: cofunc = self.dual_arg @@ -626,9 +628,9 @@ def callable() -> Function | Number: else: f.dat.data_wo[:] = f_point_eval_input_ordering.dat.data_ro[:] - if into_quadrature_space: - f_target = Function(self.target_space) - assemble(interpolate(f, self.target_space), tensor=f_target) + if self.into_quadrature_space: + f_target = Function(self.original_target_space) + assemble(interpolate(f, self.original_target_space), tensor=f_target) else: f_target = f From 35589e7cda8c73e5e62be82ba092f0cd11d41fb0 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Thu, 5 Feb 2026 13:14:14 +0000 Subject: [PATCH 26/34] lint --- firedrake/interpolation.py | 8 ++++---- .../firedrake/regression/test_cross_mesh_non_lagrange.py | 5 ++--- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 10758a9477..4f063a19a7 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -430,15 +430,15 @@ def __init__(self, expr: Interpolate): if self.source_mesh.unique().geometric_dimension != self.target_mesh.unique().geometric_dimension: raise ValueError("Geometric dimensions of source and destination meshes must match.") - + # Interpolate into intermediate quadrature space for non-identity mapped elements if into_quadrature_space := not self.target_space.finat_element.has_pointwise_dual_basis: self.original_target_space = self.target_space - r"""The original target space for interpolation, as specified by the user. + r"""The original target space for interpolation, as specified by the user. This is only used if ``self.into_quadrature_space`` is ``True``.""" self.target_space = self.target_space.quadrature_space() - r"""The target space for the cross-mesh interpolation. Must have point-evaluation dofs. - If ``self.original_target_space`` does not have point-evaluation dofs, then this is + r"""The target space for the cross-mesh interpolation. Must have point-evaluation dofs. + If ``self.original_target_space`` does not have point-evaluation dofs, then this is an intermediate quadrature space.""" self.into_quadrature_space = into_quadrature_space diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index e0daba122e..01890852bd 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -1,7 +1,6 @@ from firedrake import * import pytest import numpy as np -from finat.quadrature import QuadratureRule from functools import partial @@ -24,9 +23,9 @@ def fs_shape(V): raise ValueError("Invalid function space shape") -@pytest.fixture(params=[("RT", 1), ("RT", 2), ("BDM", 1), ("BDM", 2), ("BDFM", 2), +@pytest.fixture(params=[("RT", 1), ("RT", 2), ("BDM", 1), ("BDM", 2), ("BDFM", 2), ("HHJ", 0), ("HHJ", 2), ("N1curl", 1), ("N1curl", 2), - ("N2curl", 1), ("N2curl", 2), ("GLS", 1), ("GLS", 2), + ("N2curl", 1), ("N2curl", 2), ("GLS", 1), ("GLS", 2), ("GLS2", 2), ("Regge", 0), ("Regge", 2)], ids=lambda x: f"{x[0]}_{x[1]}") def V(request): From 48861a21ac503a440aeb45b9bd7708a742821899 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Thu, 5 Feb 2026 13:19:13 +0000 Subject: [PATCH 27/34] docstring fix --- firedrake/interpolation.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 4f063a19a7..b34483cb92 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -445,18 +445,13 @@ def __init__(self, expr: Interpolate): @cached_property def _target_space_element(self) -> FiniteElementBase: - """Return the element of `self.target_space`. If V is tensor/vector valued, - return the base scalar element. - - Parameters - ---------- - V - A :class:`.WithGeometry` function space. + """The element of `self.target_space`. If `self.target_space` is tensor/vector valued, + the base scalar element. Returns ------- FiniteElementBase - The base element of V. + The base element of `self.target_space`. """ dest_element = self.target_space.ufl_element() if isinstance(dest_element, MixedElement): @@ -490,7 +485,7 @@ def _target_space_type(self) -> Callable[..., WithGeometry]: @cached_property def _symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: - """Return symbolic ``Interpolate`` expressions for point evaluation of `self.target_space`s + """The symbolic ``Interpolate`` expressions for point evaluation of `self.target_space`s dofs in the source mesh, and the corresponding input-ordering interpolation. Returns From 6d94bb25723b79a5c19555cf4e6f487f3b65d6c3 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 9 Feb 2026 16:59:27 +0000 Subject: [PATCH 28/34] conflicts --- firedrake/interpolation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index b34483cb92..cd0332a9c7 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -19,7 +19,7 @@ from pyop2 import op2 from pyop2.caching import memory_and_disk_cache -from finat.ufl import TensorElement, VectorElement, MixedElement +from finat.ufl import TensorElement, VectorElement, MixedElement, FiniteElementBase from finat.element_factory import create_element from tsfc.driver import compile_expression_dual_evaluation From f795fe82ca1bff37cf46be64c0eb0d0acd17b258 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 9 Feb 2026 16:59:29 +0000 Subject: [PATCH 29/34] docs --- docs/source/interpolation.rst | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/docs/source/interpolation.rst b/docs/source/interpolation.rst index 116b34cbef..50441902ff 100644 --- a/docs/source/interpolation.rst +++ b/docs/source/interpolation.rst @@ -261,17 +261,11 @@ which indeed contracts into a number. Interpolation across meshes --------------------------- -The interpolation API supports interpolation between meshes where the target -function space has finite elements (as given in the list of -:ref:`supported elements `) - -* **Lagrange/CG** (also known as Continuous Galerkin or P elements), -* **Q** (i.e. Lagrange/CG on lines, quadrilaterals and hexahedra), -* **Discontinuous Lagrange/DG** (also known as Discontinuous Galerkin or DP elements) and -* **DQ** (i.e. Discontinuous Lagrange/DG on lines, quadrilaterals and hexahedra). - -Vector, tensor, and mixed function spaces can also be interpolated into from -other meshes as long as they are constructed from these spaces. +The interpolation API supports interpolation across meshes where the target +function space has any finite element which supports interpolation, as specified in the list of +:ref:`supported elements `. Vector, tensor, and mixed function +spaces can also be interpolated into from other meshes as long as they are +constructed from these spaces. .. note:: From cdeb71b1a581d17dc04fdcbfa2272c06e17aa548 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 10 Feb 2026 11:49:09 +0000 Subject: [PATCH 30/34] fixes --- firedrake/interpolation.py | 14 ++++++-------- .../regression/test_cross_mesh_non_lagrange.py | 4 ++-- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index cd0332a9c7..b65e336e7c 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -434,8 +434,7 @@ def __init__(self, expr: Interpolate): # Interpolate into intermediate quadrature space for non-identity mapped elements if into_quadrature_space := not self.target_space.finat_element.has_pointwise_dual_basis: self.original_target_space = self.target_space - r"""The original target space for interpolation, as specified by the user. - This is only used if ``self.into_quadrature_space`` is ``True``.""" + r"""The original target space for interpolation, as specified by the user.""" self.target_space = self.target_space.quadrature_space() r"""The target space for the cross-mesh interpolation. Must have point-evaluation dofs. If ``self.original_target_space`` does not have point-evaluation dofs, then this is @@ -518,21 +517,21 @@ def _symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: missing_points_behaviour=self.missing_points_behaviour, ) except VertexOnlyMeshMissingPointsError: - raise DofNotDefinedError(f"The given target function space on domain {self.target_mesh} " + raise DofNotDefinedError(f"The given target function space on domain {target_mesh} " "contains degrees of freedom which cannot cannot be defined in the " - f"source function space on domain {self.source_mesh}. " + f"source function space on domain {self.source_mesh.unique()}. " "This may be because the target mesh covers a larger domain than the " "source mesh. To disable this error, set allow_missing_dofs=True.") - # Get expression for point evaluation at the dest_node_coords + # Expression for point evaluation at the dest_node_coords P0DG_vom = self._target_space_type(vom, "DG", 0) point_eval = interpolate(self.operand, P0DG_vom) - # Interpolate into the input-ordering VOM + # Expression for interpolating into the input-ordering VOM P0DG_vom_input_ordering = self._target_space_type(vom.input_ordering, "DG", 0) - arg = Argument(P0DG_vom, 0 if self.ufl_interpolate.is_adjoint else 1) point_eval_input_ordering = interpolate(arg, P0DG_vom_input_ordering) + return point_eval, point_eval_input_ordering def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None): @@ -541,7 +540,6 @@ def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None) raise NotImplementedError("bcs not implemented for cross-mesh interpolation.") mat_type = mat_type or "aij" - # Interpolate into intermediate quadrature space for non-identity mapped elements if self.into_quadrature_space: f = Function(self.target_space.dual() if self.ufl_interpolate.is_adjoint else self.target_space) else: diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 01890852bd..1ce8d9316d 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -34,7 +34,7 @@ def V(request): return FunctionSpace(mesh, element, degree) -# @pytest.mark.parallel([1, 3]) +@pytest.mark.parallel([1, 3]) @pytest.mark.parametrize("rank", [1, 2]) def test_cross_mesh(V, rank): mesh1 = UnitSquareMesh(5, 5) @@ -94,7 +94,7 @@ def test_cross_mesh(V, rank): assert np.allclose(f_interpolated_direct.dat.data_ro, f_direct.dat.data_ro) -# @pytest.mark.parallel([1, 3]) +@pytest.mark.parallel([1, 3]) @pytest.mark.parametrize("rank", [0, 1, 2]) def test_cross_mesh_adjoint(V, rank): # Can already do Lagrange -> RT adjoint From cc07ec7399abbaaae355d1ac69b23594a415e7b6 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 10 Feb 2026 11:54:40 +0000 Subject: [PATCH 31/34] cache intermediate interpolation expressions --- firedrake/interpolation.py | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index b65e336e7c..75135e4cce 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -533,6 +533,24 @@ def _symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: point_eval_input_ordering = interpolate(arg, P0DG_vom_input_ordering) return point_eval, point_eval_input_ordering + + @cached_property + def _interpolate_from_quadrature(self) -> Interpolate: + """Returns symbolic expression for interpolation from the intermediate quadrature + space into the user-provided target space. Only relevant if `self.into_quadrature_space` is True. + + Returns + ------- + Interpolate + A symbolic interpolate expression. + """ + if self.rank == 2: + if self.ufl_interpolate.is_adjoint: + return interpolate(TestFunction(self.target_space), self.original_target_space) + else: + return interpolate(TrialFunction(self.target_space), self.original_target_space) + elif self.ufl_interpolate.is_adjoint: + return interpolate(TestFunction(self.target_space), self.dual_arg) def _get_callable(self, tensor=None, bcs=None, mat_type=None, sub_mat_type=None): from firedrake.assemble import assemble @@ -564,10 +582,10 @@ def callable() -> PETSc.Mat: source_space = self.ufl_interpolate.function_space() if self.ufl_interpolate.is_adjoint: I = AssembledMatrix((Argument(source_space, 0), Argument(self.target_space.dual(), 1)), None, res) - return assemble(action(I, interpolate(TestFunction(self.target_space), self.original_target_space))).petscmat + return assemble(action(I, self._interpolate_from_quadrature)).petscmat else: I = AssembledMatrix((Argument(self.target_space.dual(), 0), Argument(source_space, 1)), None, res) - return assemble(action(interpolate(TrialFunction(self.target_space), self.original_target_space), I)).petscmat + return assemble(action(self._interpolate_from_quadrature, I)).petscmat else: return res @@ -576,7 +594,7 @@ def callable() -> PETSc.Mat: def callable() -> Cofunction: if self.into_quadrature_space: - cofunc = assemble(interpolate(TestFunction(self.target_space), self.dual_arg)) + cofunc = assemble(self._interpolate_from_quadrature) f_target = Cofunction(point_eval.function_space()) else: cofunc = self.dual_arg From bf6cb13cf8dbc60732883b620cd8f607ba24b7bc Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 10 Feb 2026 15:29:20 +0000 Subject: [PATCH 32/34] review --- firedrake/interpolation.py | 6 ++-- .../test_cross_mesh_non_lagrange.py | 32 ++++++++++++------- 2 files changed, 25 insertions(+), 13 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 75135e4cce..a5304d7b6c 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -431,7 +431,7 @@ def __init__(self, expr: Interpolate): if self.source_mesh.unique().geometric_dimension != self.target_mesh.unique().geometric_dimension: raise ValueError("Geometric dimensions of source and destination meshes must match.") - # Interpolate into intermediate quadrature space for non-identity mapped elements + # Interpolate into intermediate quadrature space for non-point-evaluation elements if into_quadrature_space := not self.target_space.finat_element.has_pointwise_dual_basis: self.original_target_space = self.target_space r"""The original target space for interpolation, as specified by the user.""" @@ -505,6 +505,8 @@ def _symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: raise DofTypeError(f"FunctionSpace {self.target_space} must have point-evaluation dofs.") # Immerse coordinates of target space point evaluation dofs in src_mesh + # If `self.into_quadrature_space` is true, then the point evaluation dofs + # are the quadrature points of the original target space. target_mesh = self.target_space.mesh().unique() target_space_vec = VectorFunctionSpace(target_mesh, self._target_space_element) f_dest_node_coords = assemble(interpolate(target_mesh.coordinates, target_space_vec)) @@ -533,7 +535,7 @@ def _symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: point_eval_input_ordering = interpolate(arg, P0DG_vom_input_ordering) return point_eval, point_eval_input_ordering - + @cached_property def _interpolate_from_quadrature(self) -> Interpolate: """Returns symbolic expression for interpolation from the intermediate quadrature diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 1ce8d9316d..40c0e14000 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -30,7 +30,7 @@ def fs_shape(V): ids=lambda x: f"{x[0]}_{x[1]}") def V(request): element, degree = request.param - mesh = UnitSquareMesh(8, 8) + mesh = UnitSquareMesh(16, 16) return FunctionSpace(mesh, element, degree) @@ -102,12 +102,11 @@ def test_cross_mesh_adjoint(V, rank): name = V.ufl_element()._short_name deg = V.ufl_element().degree() if name in ["N1curl", "GLS", "RT"] and deg == 1: - if name == "RT" and rank == 0: - pass - else: - pytest.skip(f"Not exact for degree {deg} {name} elements") + exact = False elif name in ["Regge", "HHJ"] and deg == 0: - pytest.skip(f"Not exact for degree {deg} {name} elements") + exact = False + else: + exact = True mesh1 = UnitSquareMesh(2, 2) x1 = SpatialCoordinate(mesh1) @@ -129,6 +128,16 @@ def test_cross_mesh_adjoint(V, rank): oneform_V = inner(expr, TestFunction(V)) * dx # V^* cofunc_Vtarget_direct = assemble(inner(target_expr, TestFunction(V_target)) * dx) + if exact: + def close(x, y): + if rank == 0: + return np.isclose(x, y) + else: + return np.allclose(x, y) + else: + def close(x, y): + return np.linalg.norm(x - y)**2 < 1e-5 + Q = V.quadrature_space() if rank == 2: @@ -143,9 +152,10 @@ def test_cross_mesh_adjoint(V, rank): assert mat_equals(I_manual, I_direct) cofunc_Vtarget_manual = assemble(action(I_manual, oneform_V)) - assert np.allclose(cofunc_Vtarget_manual.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) + assert close(cofunc_Vtarget_manual.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) + cofunc_Vtarget = assemble(action(I_direct, oneform_V)) - assert np.allclose(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) + assert close(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) elif rank == 1: # Interp V^* -> Q^* I1_adj = interpolate(TestFunction(Q), oneform_V) # SameMesh @@ -154,11 +164,11 @@ def test_cross_mesh_adjoint(V, rank): # Interp Q^* -> V_target^* I2_adj = interpolate(TestFunction(V_target), cofunc_Q) # CrossMesh cofunc_Vtarget_manual = assemble(I2_adj) - assert np.allclose(cofunc_Vtarget_manual.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) + assert close(cofunc_Vtarget_manual.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) cofunc_Vtarget = assemble(interpolate(TestFunction(V_target), oneform_V)) # V^* -> V_target^* - assert np.allclose(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) + assert close(cofunc_Vtarget.dat.data_ro, cofunc_Vtarget_direct.dat.data_ro) elif rank == 0: res = assemble(interpolate(target_expr, oneform_V)) actual = assemble(inner(expr, expr) * dx) - assert np.isclose(res, actual) + assert close(res, actual) From a2cce69fb6981deeef627c47635612b1b0f33d82 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 11 Feb 2026 14:10:16 +0000 Subject: [PATCH 33/34] review --- tests/firedrake/regression/test_cross_mesh_non_lagrange.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index 40c0e14000..e1c22d9aa7 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -94,7 +94,7 @@ def test_cross_mesh(V, rank): assert np.allclose(f_interpolated_direct.dat.data_ro, f_direct.dat.data_ro) -@pytest.mark.parallel([1, 3]) +# @pytest.mark.parallel([1, 3]) @pytest.mark.parametrize("rank", [0, 1, 2]) def test_cross_mesh_adjoint(V, rank): # Can already do Lagrange -> RT adjoint @@ -136,7 +136,7 @@ def close(x, y): return np.allclose(x, y) else: def close(x, y): - return np.linalg.norm(x - y)**2 < 1e-5 + return np.linalg.norm(x - y) < 0.003 Q = V.quadrature_space() From 726e454cd72d124d2d98aeb009e9f0ce4973f6d7 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 11 Feb 2026 14:15:45 +0000 Subject: [PATCH 34/34] ............... --- tests/firedrake/regression/test_cross_mesh_non_lagrange.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py index e1c22d9aa7..ac0ea40ecb 100644 --- a/tests/firedrake/regression/test_cross_mesh_non_lagrange.py +++ b/tests/firedrake/regression/test_cross_mesh_non_lagrange.py @@ -94,7 +94,7 @@ def test_cross_mesh(V, rank): assert np.allclose(f_interpolated_direct.dat.data_ro, f_direct.dat.data_ro) -# @pytest.mark.parallel([1, 3]) +@pytest.mark.parallel([1, 3]) @pytest.mark.parametrize("rank", [0, 1, 2]) def test_cross_mesh_adjoint(V, rank): # Can already do Lagrange -> RT adjoint