Add cross mesh test

fuse/spaces/element_sobolev_spaces.py

@@ -1,4 +1,5 @@
 from functools import total_ordering
+from ufl.sobolevspace import H1, HDiv, HCurl, L2, H2
 
 
 @total_ordering
@@ -55,6 +56,9 @@ def __init__(self, cell):
     def __repr__(self):
         return "H1"
 
+    def to_ufl(self):
+        return H1
+
 
 class CellHDiv(ElementSobolevSpace):
 
@@ -64,6 +68,9 @@ def __init__(self, cell):
     def __repr__(self):
         return "HDiv"
 
+    def to_ufl(self):
+        return HDiv
+
 
 class CellHCurl(ElementSobolevSpace):
 
@@ -73,6 +80,9 @@ def __init__(self, cell):
     def __repr__(self):
         return "HCurl"
 
+    def to_ufl(self):
+        return HCurl
+
 
 class CellH2(ElementSobolevSpace):
 
@@ -82,6 +92,9 @@ def __init__(self, cell):
     def __repr__(self):
         return "H2"
 
+    def to_ufl(self):
+        return H2
+
 
 class CellL2(ElementSobolevSpace):
 
@@ -90,3 +103,6 @@ def __init__(self, cell):
 
     def __repr__(self):
         return "L2"
+
+    def to_ufl(self):
+        return L2

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,25 +1073,27 @@ 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),
                                                             (construct_tet_ned_2nd_kind_2, "N2curl", 2, 1e-12),
                                                             (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)
-    ufl_elem1 = elem1.to_ufl()
+    if hasattr(elem_gen, "__call__"):
+        elem = elem_gen()
+    else:
+        elem = elem_gen
+    ufl_elem = elem.to_ufl()
 
     def expr(mesh):
         x = SpatialCoordinate(mesh)
@@ -1098,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)

test/test_interpolation.py

@@ -0,0 +1,92 @@
+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_tri_to_quad():
+    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)))
+
+
+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)