From 11dc1338a92d608764d75b5011d2607cb4bb489d Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 10 Dec 2024 17:42:28 +0000 Subject: [PATCH] Add test/finat/conftest.py (#114) * Add conftest.py * set --import-mode=importlib * Extend MyMapping to 3D --- pyproject.toml | 3 + test/finat/conftest.py | 128 +++++++++++++++++++++++++++ test/finat/fiat_mapping.py | 72 --------------- test/finat/test_mass_conditioning.py | 14 +-- test/finat/test_zany_mapping.py | 54 ++++------- 5 files changed, 152 insertions(+), 119 deletions(-) create mode 100644 test/finat/conftest.py delete mode 100644 test/finat/fiat_mapping.py diff --git a/pyproject.toml b/pyproject.toml index 6e10c7b8a..c224c1fbe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -37,3 +37,6 @@ test = ["pytest"] [tool.setuptools] packages = ["FIAT", "finat", "finat.ufl", "gem"] + +[tool.pytest.ini_options] +addopts = "--import-mode=importlib" diff --git a/test/finat/conftest.py b/test/finat/conftest.py new file mode 100644 index 000000000..c31ef7350 --- /dev/null +++ b/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} diff --git a/test/finat/fiat_mapping.py b/test/finat/fiat_mapping.py deleted file mode 100644 index d857b2577..000000000 --- a/test/finat/fiat_mapping.py +++ /dev/null @@ -1,72 +0,0 @@ -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): - # Firedrake interprets this as 2x the circumradius - # cs = (np.prod([self.phys_cell.volume_of_subcomplex(1, i) - # for i in range(3)]) - # / 2.0 / self.phys_cell.volume()) - # return np.asarray([cs for _ in range(3)]) - # Currently, just return 1 so we can compare FIAT dofs - # to transformed dofs. - - return np.ones((3,)) - - 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): - return gem.Literal(np.asarray([self.ref_cell.compute_normalized_edge_tangent(i) for i in range(3)])) - - def reference_normals(self): - return gem.Literal( - np.asarray([self.ref_cell.compute_normal(i) - for i in range(3)])) - - def physical_normals(self): - return gem.Literal( - np.asarray([self.phys_cell.compute_normal(i) - for i in range(3)])) - - def physical_tangents(self): - return gem.Literal( - np.asarray([self.phys_cell.compute_normalized_edge_tangent(i) - for i in range(3)])) - - def physical_edge_lengths(self): - return gem.Literal( - np.asarray([self.phys_cell.volume_of_subcomplex(1, i) - for i in range(3)])) - - 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 FiredrakeMapping(MyMapping): - - def cell_size(self): - # Firedrake interprets this as 2x the circumradius - cs = (np.prod([self.phys_cell.volume_of_subcomplex(1, i) - for i in range(3)]) - / 2.0 / self.phys_cell.volume()) - return np.asarray([cs for _ in range(3)]) diff --git a/test/finat/test_mass_conditioning.py b/test/finat/test_mass_conditioning.py index cc4eafcaa..c6de89cc9 100644 --- a/test/finat/test_mass_conditioning.py +++ b/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 diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 9444f5e24..d5136e651 100644 --- a/test/finat/test_zany_mapping.py +++ b/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,30 +85,13 @@ 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", [ @@ -117,23 +99,23 @@ def test_C1_elements(ref_el, phys_el, element): 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)