Add test/finat/conftest.py (#114)

* Add conftest.py

* set --import-mode=importlib

* Extend MyMapping to 3D

pyproject.toml

@@ -37,3 +37,6 @@ test = ["pytest"]
 
 [tool.setuptools]
 packages = ["FIAT", "finat", "finat.ufl", "gem"]
+
+[tool.pytest.ini_options]
+addopts = "--import-mode=importlib"

test/finat/conftest.py

@@ -0,0 +1,128 @@
+import pytest
+import FIAT
+import gem
+import numpy as np
+from finat.physically_mapped import PhysicalGeometry
+
+
+class MyMapping(PhysicalGeometry):
+    def __init__(self, ref_cell, phys_cell):
+        self.ref_cell = ref_cell
+        self.phys_cell = phys_cell
+
+        self.A, self.b = FIAT.reference_element.make_affine_mapping(
+            self.ref_cell.vertices,
+            self.phys_cell.vertices)
+
+    def cell_size(self):
+        # Currently, just return 1 so we can compare FIAT dofs
+        # to transformed dofs.
+        return np.ones((len(self.ref_cell.vertices),))
+
+    def detJ_at(self, point):
+        return gem.Literal(np.linalg.det(self.A))
+
+    def jacobian_at(self, point):
+        return gem.Literal(self.A)
+
+    def normalized_reference_edge_tangents(self):
+        top = self.ref_cell.get_topology()
+        return gem.Literal(
+            np.asarray([self.ref_cell.compute_normalized_edge_tangent(i)
+                        for i in sorted(top[1])]))
+
+    def reference_normals(self):
+        sd = self.ref_cell.get_spatial_dimension()
+        top = self.ref_cell.get_topology()
+        return gem.Literal(
+            np.asarray([self.ref_cell.compute_normal(i)
+                        for i in sorted(top[sd-1])]))
+
+    def physical_normals(self):
+        sd = self.phys_cell.get_spatial_dimension()
+        top = self.phys_cell.get_topology()
+        return gem.Literal(
+            np.asarray([self.phys_cell.compute_normal(i)
+                        for i in sorted(top[sd-1])]))
+
+    def physical_tangents(self):
+        top = self.phys_cell.get_topology()
+        return gem.Literal(
+            np.asarray([self.phys_cell.compute_normalized_edge_tangent(i)
+                        for i in sorted(top[1])]))
+
+    def physical_edge_lengths(self):
+        top = self.phys_cell.get_topology()
+        return gem.Literal(
+            np.asarray([self.phys_cell.volume_of_subcomplex(1, i)
+                        for i in sorted(top[1])]))
+
+    def physical_points(self, ps, entity=None):
+        prefs = ps.points
+        A, b = self.A, self.b
+        return gem.Literal(np.asarray([A @ x + b for x in prefs]))
+
+    def physical_vertices(self):
+        return gem.Literal(self.phys_cell.verts)
+
+
+class ScaledMapping(MyMapping):
+
+    def cell_size(self):
+        # Firedrake interprets this as 2x the circumradius
+        sd = self.phys_cell.get_spatial_dimension()
+        top = self.phys_cell.get_topology()
+        vol = self.phys_cell.volume()
+        edges = [self.phys_cell.volume_of_subcomplex(1, i)
+                 for i in sorted(top[1])]
+
+        if sd == 1:
+            cs = vol
+        elif sd == 2:
+            cs = np.prod(edges) / (2 * vol)
+        elif sd == 3:
+            edge_pairs = [edges[i] * edges[j] for i in top[1] for j in top[1]
+                          if len(set(top[1][i] + top[1][j])) == len(top[0])]
+            cs = 1.0 / (12 * vol)
+            for k in range(4):
+                s = [1]*len(edge_pairs)
+                if k > 0:
+                    s[k-1] = -1
+                cs *= np.dot(s, edge_pairs) ** 0.5
+        else:
+            raise NotImplementedError(f"Cell size not implemented in {sd} dimensions")
+        return np.asarray([cs for _ in sorted(top[0])])
+
+
+def scaled_simplex(dim, scale):
+    K = FIAT.ufc_simplex(dim)
+    K.vertices = scale * np.array(K.vertices)
+    return K
+
+
+@pytest.fixture
+def ref_el():
+    K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
+    return K
+
+
+@pytest.fixture
+def phys_el():
+    K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
+    K[2].vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84))
+    K[3].vertices = ((0, 0, 0),
+                     (1., 0.1, -0.37),
+                     (0.01, 0.987, -.23),
+                     (-0.1, -0.2, 1.38))
+    return K
+
+
+@pytest.fixture
+def ref_to_phys(ref_el, phys_el):
+    return {dim: MyMapping(ref_el[dim], phys_el[dim]) for dim in ref_el}
+
+
+@pytest.fixture
+def scaled_ref_to_phys(ref_el):
+    return {dim: [ScaledMapping(ref_el[dim], scaled_simplex(dim, 0.5**k)) for k in range(3)]
+            for dim in ref_el}

test/finat/test_mass_conditioning.py

@@ -1,11 +1,8 @@
-import FIAT
 import finat
 import numpy as np
 import pytest
 from gem.interpreter import evaluate
 
-from fiat_mapping import FiredrakeMapping
-
 
 @pytest.mark.parametrize("element, degree, variant", [
     (finat.Hermite, 3, None),
@@ -19,9 +16,9 @@
     (finat.Argyris, 5, None),
     (finat.Argyris, 6, None),
 ])
-def test_mass_scaling(element, degree, variant):
+def test_mass_scaling(scaled_ref_to_phys, element, degree, variant):
     sd = 2
-    ref_cell = FIAT.ufc_simplex(sd)
+    ref_cell = scaled_ref_to_phys[sd][0].ref_cell
     if variant is not None:
         ref_element = element(ref_cell, degree, variant=variant)
     else:
@@ -32,12 +29,7 @@ def test_mass_scaling(element, degree, variant):
     qwts = Q.weights
 
     kappa = []
-    for k in range(3):
-        h = 2 ** -k
-        phys_cell = FIAT.ufc_simplex(2)
-        new_verts = h * np.array(phys_cell.get_vertices())
-        phys_cell.vertices = tuple(map(tuple, new_verts))
-        mapping = FiredrakeMapping(ref_cell, phys_cell)
+    for mapping in scaled_ref_to_phys[sd]:
         J_gem = mapping.jacobian_at(ref_cell.make_points(sd, 0, sd+1)[0])
         J = evaluate([J_gem])[0].arr
 

test/finat/test_zany_mapping.py

@@ -4,8 +4,6 @@
 import pytest
 from gem.interpreter import evaluate
 
-from fiat_mapping import MyMapping
-
 
 def make_unisolvent_points(element, interior=False):
     degree = element.degree()
@@ -23,7 +21,9 @@ def make_unisolvent_points(element, interior=False):
     return pts
 
 
-def check_zany_mapping(element, ref_cell, phys_cell, *args, **kwargs):
+def check_zany_mapping(element, ref_to_phys, *args, **kwargs):
+    phys_cell = ref_to_phys.phys_cell
+    ref_cell = ref_to_phys.ref_cell
     phys_element = element(phys_cell, *args, **kwargs).fiat_equivalent
     finat_element = element(ref_cell, *args, **kwargs)
 
@@ -65,9 +65,8 @@ def check_zany_mapping(element, ref_cell, phys_cell, *args, **kwargs):
     # Zany map the results
     num_bfs = phys_element.space_dimension()
     num_dofs = finat_element.space_dimension()
-    mappng = MyMapping(ref_cell, phys_cell)
     try:
-        Mgem = finat_element.basis_transformation(mappng)
+        Mgem = finat_element.basis_transformation(ref_to_phys)
         M = evaluate([Mgem])[0].arr
         ref_vals_zany = np.tensordot(M, ref_vals_piola, (-1, 0))
     except AttributeError:
@@ -86,54 +85,37 @@ def check_zany_mapping(element, ref_cell, phys_cell, *args, **kwargs):
     assert np.allclose(ref_vals_zany, phys_vals[:num_dofs])
 
 
-@pytest.fixture
-def ref_el(request):
-    K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
-    return K
-
-
-@pytest.fixture
-def phys_el(request):
-    K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
-    K[2].vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84))
-    K[3].vertices = ((0, 0, 0),
-                     (1., 0.1, -0.37),
-                     (0.01, 0.987, -.23),
-                     (-0.1, -0.2, 1.38))
-    return K
-
-
 @pytest.mark.parametrize("element", [
                          finat.Morley,
                          finat.Hermite,
                          finat.Bell,
                          ])
-def test_C1_elements(ref_el, phys_el, element):
-    check_zany_mapping(element, ref_el[2], phys_el[2])
+def test_C1_elements(ref_to_phys, element):
+    check_zany_mapping(element, ref_to_phys[2])
 
 
 @pytest.mark.parametrize("element", [
                          finat.QuadraticPowellSabin6,
                          finat.QuadraticPowellSabin12,
                          finat.ReducedHsiehCloughTocher,
                          ])
-def test_C1_macroelements(ref_el, phys_el, element):
+def test_C1_macroelements(ref_to_phys, element):
     kwargs = {}
     if element == finat.QuadraticPowellSabin12:
         kwargs = dict(avg=True)
-    check_zany_mapping(element, ref_el[2], phys_el[2], **kwargs)
+    check_zany_mapping(element, ref_to_phys[2], **kwargs)
 
 
 @pytest.mark.parametrize("element, degree", [
     *((finat.Argyris, k) for k in range(5, 8)),
     *((finat.HsiehCloughTocher, k) for k in range(3, 6))
 ])
-def test_high_order_C1_elements(ref_el, phys_el, element, degree):
-    check_zany_mapping(element, ref_el[2], phys_el[2], degree, avg=True)
+def test_high_order_C1_elements(ref_to_phys, element, degree):
+    check_zany_mapping(element, ref_to_phys[2], degree, avg=True)
 
 
-def test_argyris_point(ref_el, phys_el):
-    check_zany_mapping(finat.Argyris, ref_el[2], phys_el[2], variant="point")
+def test_argyris_point(ref_to_phys):
+    check_zany_mapping(finat.Argyris, ref_to_phys[2], variant="point")
 
 
 zany_piola_elements = {
@@ -162,15 +144,15 @@ def test_argyris_point(ref_el, phys_el):
                          *((2, e) for e in zany_piola_elements[3]),
                          *((3, e) for e in zany_piola_elements[3]),
                          ])
-def test_piola(ref_el, phys_el, element, dimension):
-    check_zany_mapping(element, ref_el[dimension], phys_el[dimension])
+def test_piola(ref_to_phys, element, dimension):
+    check_zany_mapping(element, ref_to_phys[dimension])
 
 
 @pytest.mark.parametrize("element, degree, variant", [
     *((finat.HuZhang, k, v) for v in ("integral", "point") for k in range(3, 6)),
 ])
-def test_piola_triangle_high_order(ref_el, phys_el, element, degree, variant):
-    check_zany_mapping(element, ref_el[2], phys_el[2], degree, variant)
+def test_piola_triangle_high_order(ref_to_phys, element, degree, variant):
+    check_zany_mapping(element, ref_to_phys[2], degree, variant)
 
 
 @pytest.mark.parametrize("element, degree", [
@@ -180,5 +162,5 @@ def test_piola_triangle_high_order(ref_el, phys_el, element, degree, variant):
                          *((finat.GopalakrishnanLedererSchoberlSecondKind, k) for k in range(0, 3)),
                          ])
 @pytest.mark.parametrize("dimension", [2, 3])
-def test_affine(ref_el, phys_el, element, degree, dimension):
-    check_zany_mapping(element, ref_el[dimension], phys_el[dimension], degree)
+def test_affine(ref_to_phys, element, degree, dimension):
+    check_zany_mapping(element, ref_to_phys[dimension], degree)