From 187558662156f9002a4aaa16846c35c30bfed6f0 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 11 Mar 2026 15:36:42 +0000 Subject: [PATCH 1/6] use ufl sobolev spaces --- fuse/spaces/element_sobolev_spaces.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/fuse/spaces/element_sobolev_spaces.py b/fuse/spaces/element_sobolev_spaces.py index 2a3aa4fb..07c9fdd9 100644 --- a/fuse/spaces/element_sobolev_spaces.py +++ b/fuse/spaces/element_sobolev_spaces.py @@ -1,5 +1,5 @@ from functools import total_ordering - +from ufl.sobolevspace import H1, HDiv, HCurl, L2, H2 @total_ordering class ElementSobolevSpace(object): @@ -54,6 +54,9 @@ def __init__(self, cell): def __repr__(self): return "H1" + + def to_ufl(self): + return H1 class CellHDiv(ElementSobolevSpace): @@ -64,6 +67,9 @@ def __init__(self, cell): def __repr__(self): return "HDiv" + def to_ufl(self): + return HDiv + class CellHCurl(ElementSobolevSpace): @@ -73,6 +79,9 @@ def __init__(self, cell): def __repr__(self): return "HCurl" + def to_ufl(self): + return HCurl + class CellH2(ElementSobolevSpace): @@ -82,6 +91,9 @@ def __init__(self, cell): def __repr__(self): return "H2" + def to_ufl(self): + return H2 + class CellL2(ElementSobolevSpace): @@ -90,3 +102,6 @@ def __init__(self, cell): def __repr__(self): return "L2" + + def to_ufl(self): + return L2 From cb43446fe1cb90cbfa6fc1c50bc4838b0426638f Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 9 Jul 2026 15:55:24 +0100 Subject: [PATCH 2/6] lint --- fuse/spaces/element_sobolev_spaces.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/fuse/spaces/element_sobolev_spaces.py b/fuse/spaces/element_sobolev_spaces.py index 07c9fdd9..b5964d9d 100644 --- a/fuse/spaces/element_sobolev_spaces.py +++ b/fuse/spaces/element_sobolev_spaces.py @@ -1,6 +1,7 @@ from functools import total_ordering from ufl.sobolevspace import H1, HDiv, HCurl, L2, H2 + @total_ordering class ElementSobolevSpace(object): """ @@ -54,7 +55,7 @@ def __init__(self, cell): def __repr__(self): return "H1" - + def to_ufl(self): return H1 From 4a0191e936d22c92f6084c67878cdebbe9e171da Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 10 Jul 2026 14:32:28 +0100 Subject: [PATCH 3/6] lint more --- test/test_convert_to_fiat.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index b90440ad..4c1f53c5 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -1076,12 +1076,14 @@ def vec(mesh): (construct_tet_ned_2nd_kind_2_non_bary, "N2curl", 2, 1e-12), (construct_tet_ned_2nd_kind_3, "N2curl", 3, 1e-12), (construct_tet_ned2, "N1curl", 2, 1e-13), - (lambda cell: periodic_table(1, 3, 1, 3), "N2curl", 3, 1e-12), - (lambda cell: periodic_table(1, 3, 1, 4), "N2curl", 4, 1e-12), + (periodic_table(1, 3, 1, 3), "N2curl", 3, 1e-12), + (periodic_table(1, 3, 1, 4), "N2curl", 4, 1e-12), (construct_tet_ned3_old, "N1curl", 2, 1e-13)]) def test_two_tet_projection(elem_gen, elem_code, deg, max_err): - cell = make_tetrahedron() - elem1 = elem_gen(cell) + if hasattr(elem_gen, "__call__"): + elem1 = elem_gen() + else: + elem1 = elem_gen ufl_elem1 = elem1.to_ufl() def expr(mesh): From a088c94966b14a08919f095e1199aabbad1f18ea Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 10 Jul 2026 16:19:54 +0100 Subject: [PATCH 4/6] add test for cross mesh interpolation --- test/test_interpolation.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 test/test_interpolation.py diff --git a/test/test_interpolation.py b/test/test_interpolation.py new file mode 100644 index 00000000..af894d84 --- /dev/null +++ b/test/test_interpolation.py @@ -0,0 +1,21 @@ +from firedrake import * +from fuse import * +import numpy as np +from test_2d_examples_docs import construct_cg1, construct_cg3 +from test_convert_to_fiat import create_cg2, create_cg2_tri + +def test_cross_mesh(): + mesh1 = UnitSquareMesh(10, 10, use_fuse=True) + mesh2 = UnitSquareMesh(10, 10, quadrilateral=True, use_fuse=True) + A = create_cg2() + B = create_cg2() + V1 = FunctionSpace(mesh1, create_cg2_tri().to_ufl()) + V2 = FunctionSpace(mesh2, tensor_product(A, B).flatten().to_ufl()) + + f1 = Function(V1) + f2 = Function(V2) + + x = SpatialCoordinate(mesh1) + f1 = f1.interpolate(x[0]**2 + x[1]**2) + f2 = f2.interpolate(f1) + assert np.allclose(sqrt(assemble(inner(f1, f1) * dx)), sqrt(assemble(inner(f2, f2) * dx))) \ No newline at end of file From 8640f8ddaba0768be49ca153ee2595765b07f07f Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 10 Jul 2026 16:23:49 +0100 Subject: [PATCH 5/6] lint --- test/test_interpolation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index af894d84..457eec19 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -1,9 +1,9 @@ from firedrake import * from fuse import * import numpy as np -from test_2d_examples_docs import construct_cg1, construct_cg3 from test_convert_to_fiat import create_cg2, create_cg2_tri + def test_cross_mesh(): mesh1 = UnitSquareMesh(10, 10, use_fuse=True) mesh2 = UnitSquareMesh(10, 10, quadrilateral=True, use_fuse=True) @@ -18,4 +18,4 @@ def test_cross_mesh(): x = SpatialCoordinate(mesh1) f1 = f1.interpolate(x[0]**2 + x[1]**2) f2 = f2.interpolate(f1) - assert np.allclose(sqrt(assemble(inner(f1, f1) * dx)), sqrt(assemble(inner(f2, f2) * dx))) \ No newline at end of file + assert np.allclose(sqrt(assemble(inner(f1, f1) * dx)), sqrt(assemble(inner(f2, f2) * dx))) From ff8ac8bf017b70f29b2f9a10be08cf0d75d2b6d0 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 13 Jul 2026 10:34:21 +0100 Subject: [PATCH 6/6] new cross mesh test and fixing others --- test/test_convert_to_fiat.py | 35 ++++++++++------- test/test_interpolation.py | 73 +++++++++++++++++++++++++++++++++++- 2 files changed, 94 insertions(+), 14 deletions(-) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 4c1f53c5..33cec066 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -147,8 +147,10 @@ def create_cg1_flipped(cell): return cg -def create_cg2(cell): +def create_cg2(cell=None): deg = 2 + if cell is None: + cell = Point(1, [Point(0), Point(0)], vertex_num=2) if cell.dim() > 1: raise NotImplementedError("This method is for cg2 on edges, please use create_cg2_tri for triangles") vert_dg = create_dg0(cell.vertices()[0]) @@ -161,8 +163,10 @@ def create_cg2(cell): return cg -def create_cg2_tri(cell): +def create_cg2_tri(cell=None): deg = 2 + if cell is None: + cell = polygon(3) Pk = PolynomialSpace(deg) vert_dg0 = create_dg0(cell.vertices()[0]) @@ -218,7 +222,9 @@ def create_cg2_tet(cell): return cg2 -def create_cg3_tet(cell, perm=True): +def create_cg3_tet(cell=None, perm=True): + if cell is None: + cell = make_tetrahedron() vert = cell.vertices()[0] edge = cell.edges()[0] @@ -1017,11 +1023,14 @@ def evaluate_pt_dict(pt_dict, fn): (construct_tet_cg4, "CG", 4), (construct_tet_cg6, "CG", 6), (construct_tet_ned3_old, "N1curl", 3), - (lambda cell: periodic_table(1, 3, 1, 3), "N2curl", 3), - (lambda cell: periodic_table(0, 3, 1, 3), "N1curl", 3)]) + (periodic_table(1, 3, 1, 3), "N2curl", 3), + (periodic_table(0, 3, 1, 3), "N1curl", 3)]) def test_two_tet_interpolation(elem_gen, elem_code, deg): cell = make_tetrahedron() - elem = elem_gen(cell) + if hasattr(elem_gen, "__call__"): + elem = elem_gen(cell) + else: + elem = elem_gen def vec(mesh): x = SpatialCoordinate(mesh) @@ -1064,11 +1073,11 @@ def vec(mesh): @pytest.mark.parametrize("elem_gen,elem_code,deg,max_err", [(construct_tet_cg6, "CG", 6, 1e-13), - (lambda cell: periodic_table(0, 3, 1, 3), "N1curl", 3, 1e-12), + (periodic_table(0, 3, 1, 3), "N1curl", 3, 1e-12), (create_cg3_tet, "CG", 3, 1e-13), (construct_tet_cg4, "CG", 4, 1e-13), - (lambda cell: periodic_table(0, 3, 0, 4), "CG", 4, 1e-13), - (lambda cell: periodic_table(0, 3, 0, 6), "CG", 6, 1e-13), + (periodic_table(0, 3, 0, 4), "CG", 4, 1e-13), + (periodic_table(0, 3, 0, 6), "CG", 6, 1e-13), (construct_tet_rt2, "RT", 2, 1e-13), (construct_tet_rt3, "RT", 3, 1e-13), (construct_tet_bdm2, "BDM", 2, 1e-13), @@ -1081,10 +1090,10 @@ def vec(mesh): (construct_tet_ned3_old, "N1curl", 2, 1e-13)]) def test_two_tet_projection(elem_gen, elem_code, deg, max_err): if hasattr(elem_gen, "__call__"): - elem1 = elem_gen() + elem = elem_gen() else: - elem1 = elem_gen - ufl_elem1 = elem1.to_ufl() + elem = elem_gen + ufl_elem = elem.to_ufl() def expr(mesh): x = SpatialCoordinate(mesh) @@ -1100,7 +1109,7 @@ def expr(mesh): sp.combinatorics.Permutation([0, 3, 2, 1]), sp.combinatorics.Permutation([0, 2, 1, 3])] - for elem in [ufl_elem1]: + for elem in [ufl_elem]: for g in group: mesh = TwoTetMesh(perm=g, use_fuse=True) print(g) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index 457eec19..7d257828 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -1,10 +1,12 @@ from firedrake import * +from firedrake.ufl_expr import extract_unique_domain from fuse import * import numpy as np +from test_2d_examples_docs import construct_cg3 from test_convert_to_fiat import create_cg2, create_cg2_tri -def test_cross_mesh(): +def test_cross_mesh_tri_to_quad(): mesh1 = UnitSquareMesh(10, 10, use_fuse=True) mesh2 = UnitSquareMesh(10, 10, quadrilateral=True, use_fuse=True) A = create_cg2() @@ -19,3 +21,72 @@ def test_cross_mesh(): f1 = f1.interpolate(x[0]**2 + x[1]**2) f2 = f2.interpolate(f1) assert np.allclose(sqrt(assemble(inner(f1, f1) * dx)), sqrt(assemble(inner(f2, f2) * dx))) + + +def test_cross_mesh_fuse_to_ufc(): + mesh1 = UnitSquareMesh(10, 10, use_fuse=True) + mesh2 = UnitSquareMesh(10, 10) + V1 = FunctionSpace(mesh1, create_cg2_tri().to_ufl()) + V2 = FunctionSpace(mesh2, "CG", 2) + + f1 = Function(V1) + f2 = Function(V2) + + x = SpatialCoordinate(mesh1) + f1 = f1.interpolate(x[0]**2 + x[1]**2) + f2 = f2.interpolate(f1) + assert np.allclose(sqrt(assemble(inner(f1, f1) * dx)), sqrt(assemble(inner(f2, f2) * dx))) + + +def test_cross_mesh(): + dest_quad = False + atol = 1e-8 + m_src = UnitSquareMesh(2, 3, use_fuse=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 + # add the coordinates of the mesh vertices to test boundaries + vertices_src = m_src.coordinates.dat.data_ro + coords = np.concatenate((coords, vertices_src)) + vertices_dest = m_dest.coordinates.dat.data_ro + coords = np.concatenate((coords, vertices_dest)) + expr_src = product(SpatialCoordinate(m_src)) + expr_dest = product(SpatialCoordinate(m_dest)) + dest_eval = PointEvaluator(m_dest, coords) + expected = np.prod(coords, axis=-1) + + V_src = FunctionSpace(m_src, construct_cg3().to_ufl()) + V_dest = FunctionSpace(m_dest, "CG", 4) + + # test_expression from test_interpolate_cross_mesh in firedrake + f_dest = assemble(interpolate(expr_src, V_dest)) + assert extract_unique_domain(f_dest) is m_dest + got = dest_eval.evaluate(f_dest) + assert np.allclose(got, expected, atol=atol) + f_dest_2 = Function(V_dest).interpolate(expr_dest) + assert np.allclose(f_dest.dat.data_ro, f_dest_2.dat.data_ro, atol=atol) + + # test_function from test_interpolate_cross_mesh in firedrake + f_dest = Function(V_dest).interpolate(expr_src) + assert extract_unique_domain(f_dest) is m_dest + + got = dest_eval.evaluate(f_dest) + assert np.allclose(got, expected, atol=atol) + + f_src = Function(V_src).interpolate(expr_src) + f_dest = assemble(interpolate(f_src, V_dest)) + assert extract_unique_domain(f_dest) is m_dest + got = dest_eval.evaluate(f_dest) + assert np.allclose(got, expected, atol=atol) + + f_dest_2 = Function(V_dest).interpolate(expr_dest) + assert np.allclose(f_dest.dat.data_ro, f_dest_2.dat.data_ro, atol=atol) + + # test Function.interpolate(...) + f_dest = Function(V_dest) + f_dest.interpolate(f_src) + assert extract_unique_domain(f_dest) is m_dest + got = dest_eval.evaluate(f_dest) + assert np.allclose(got, expected, atol=atol) + assert np.allclose(f_dest.dat.data_ro, f_dest_2.dat.data_ro, atol=atol)