Skip to content

Commit 11dc133

Browse files
pbrubeckpbrubeck
authored
Add test/finat/conftest.py (#114)
* Add conftest.py * set --import-mode=importlib * Extend MyMapping to 3D
1 parent 45f6d9e commit 11dc133

5 files changed

+152
-119
lines changed

‎pyproject.toml

+3
Original file line numberDiff line numberDiff line change
@@ -37,3 +37,6 @@ test = ["pytest"]
3737

3838
[tool.setuptools]
3939
packages = ["FIAT", "finat", "finat.ufl", "gem"]
40+
41+
[tool.pytest.ini_options]
42+
addopts = "--import-mode=importlib"

‎test/finat/conftest.py

+128
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,128 @@
1+
import pytest
2+
import FIAT
3+
import gem
4+
import numpy as np
5+
from finat.physically_mapped import PhysicalGeometry
6+
7+
8+
class MyMapping(PhysicalGeometry):
9+
def __init__(self, ref_cell, phys_cell):
10+
self.ref_cell = ref_cell
11+
self.phys_cell = phys_cell
12+
13+
self.A, self.b = FIAT.reference_element.make_affine_mapping(
14+
self.ref_cell.vertices,
15+
self.phys_cell.vertices)
16+
17+
def cell_size(self):
18+
# Currently, just return 1 so we can compare FIAT dofs
19+
# to transformed dofs.
20+
return np.ones((len(self.ref_cell.vertices),))
21+
22+
def detJ_at(self, point):
23+
return gem.Literal(np.linalg.det(self.A))
24+
25+
def jacobian_at(self, point):
26+
return gem.Literal(self.A)
27+
28+
def normalized_reference_edge_tangents(self):
29+
top = self.ref_cell.get_topology()
30+
return gem.Literal(
31+
np.asarray([self.ref_cell.compute_normalized_edge_tangent(i)
32+
for i in sorted(top[1])]))
33+
34+
def reference_normals(self):
35+
sd = self.ref_cell.get_spatial_dimension()
36+
top = self.ref_cell.get_topology()
37+
return gem.Literal(
38+
np.asarray([self.ref_cell.compute_normal(i)
39+
for i in sorted(top[sd-1])]))
40+
41+
def physical_normals(self):
42+
sd = self.phys_cell.get_spatial_dimension()
43+
top = self.phys_cell.get_topology()
44+
return gem.Literal(
45+
np.asarray([self.phys_cell.compute_normal(i)
46+
for i in sorted(top[sd-1])]))
47+
48+
def physical_tangents(self):
49+
top = self.phys_cell.get_topology()
50+
return gem.Literal(
51+
np.asarray([self.phys_cell.compute_normalized_edge_tangent(i)
52+
for i in sorted(top[1])]))
53+
54+
def physical_edge_lengths(self):
55+
top = self.phys_cell.get_topology()
56+
return gem.Literal(
57+
np.asarray([self.phys_cell.volume_of_subcomplex(1, i)
58+
for i in sorted(top[1])]))
59+
60+
def physical_points(self, ps, entity=None):
61+
prefs = ps.points
62+
A, b = self.A, self.b
63+
return gem.Literal(np.asarray([A @ x + b for x in prefs]))
64+
65+
def physical_vertices(self):
66+
return gem.Literal(self.phys_cell.verts)
67+
68+
69+
class ScaledMapping(MyMapping):
70+
71+
def cell_size(self):
72+
# Firedrake interprets this as 2x the circumradius
73+
sd = self.phys_cell.get_spatial_dimension()
74+
top = self.phys_cell.get_topology()
75+
vol = self.phys_cell.volume()
76+
edges = [self.phys_cell.volume_of_subcomplex(1, i)
77+
for i in sorted(top[1])]
78+
79+
if sd == 1:
80+
cs = vol
81+
elif sd == 2:
82+
cs = np.prod(edges) / (2 * vol)
83+
elif sd == 3:
84+
edge_pairs = [edges[i] * edges[j] for i in top[1] for j in top[1]
85+
if len(set(top[1][i] + top[1][j])) == len(top[0])]
86+
cs = 1.0 / (12 * vol)
87+
for k in range(4):
88+
s = [1]*len(edge_pairs)
89+
if k > 0:
90+
s[k-1] = -1
91+
cs *= np.dot(s, edge_pairs) ** 0.5
92+
else:
93+
raise NotImplementedError(f"Cell size not implemented in {sd} dimensions")
94+
return np.asarray([cs for _ in sorted(top[0])])
95+
96+
97+
def scaled_simplex(dim, scale):
98+
K = FIAT.ufc_simplex(dim)
99+
K.vertices = scale * np.array(K.vertices)
100+
return K
101+
102+
103+
@pytest.fixture
104+
def ref_el():
105+
K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
106+
return K
107+
108+
109+
@pytest.fixture
110+
def phys_el():
111+
K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
112+
K[2].vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84))
113+
K[3].vertices = ((0, 0, 0),
114+
(1., 0.1, -0.37),
115+
(0.01, 0.987, -.23),
116+
(-0.1, -0.2, 1.38))
117+
return K
118+
119+
120+
@pytest.fixture
121+
def ref_to_phys(ref_el, phys_el):
122+
return {dim: MyMapping(ref_el[dim], phys_el[dim]) for dim in ref_el}
123+
124+
125+
@pytest.fixture
126+
def scaled_ref_to_phys(ref_el):
127+
return {dim: [ScaledMapping(ref_el[dim], scaled_simplex(dim, 0.5**k)) for k in range(3)]
128+
for dim in ref_el}

‎test/finat/fiat_mapping.py

-72
This file was deleted.

‎test/finat/test_mass_conditioning.py

+3-11
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,8 @@
1-
import FIAT
21
import finat
32
import numpy as np
43
import pytest
54
from gem.interpreter import evaluate
65

7-
from fiat_mapping import FiredrakeMapping
8-
96

107
@pytest.mark.parametrize("element, degree, variant", [
118
(finat.Hermite, 3, None),
@@ -19,9 +16,9 @@
1916
(finat.Argyris, 5, None),
2017
(finat.Argyris, 6, None),
2118
])
22-
def test_mass_scaling(element, degree, variant):
19+
def test_mass_scaling(scaled_ref_to_phys, element, degree, variant):
2320
sd = 2
24-
ref_cell = FIAT.ufc_simplex(sd)
21+
ref_cell = scaled_ref_to_phys[sd][0].ref_cell
2522
if variant is not None:
2623
ref_element = element(ref_cell, degree, variant=variant)
2724
else:
@@ -32,12 +29,7 @@ def test_mass_scaling(element, degree, variant):
3229
qwts = Q.weights
3330

3431
kappa = []
35-
for k in range(3):
36-
h = 2 ** -k
37-
phys_cell = FIAT.ufc_simplex(2)
38-
new_verts = h * np.array(phys_cell.get_vertices())
39-
phys_cell.vertices = tuple(map(tuple, new_verts))
40-
mapping = FiredrakeMapping(ref_cell, phys_cell)
32+
for mapping in scaled_ref_to_phys[sd]:
4133
J_gem = mapping.jacobian_at(ref_cell.make_points(sd, 0, sd+1)[0])
4234
J = evaluate([J_gem])[0].arr
4335

‎test/finat/test_zany_mapping.py

+18-36
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,6 @@
44
import pytest
55
from gem.interpreter import evaluate
66

7-
from fiat_mapping import MyMapping
8-
97

108
def make_unisolvent_points(element, interior=False):
119
degree = element.degree()
@@ -23,7 +21,9 @@ def make_unisolvent_points(element, interior=False):
2321
return pts
2422

2523

26-
def check_zany_mapping(element, ref_cell, phys_cell, *args, **kwargs):
24+
def check_zany_mapping(element, ref_to_phys, *args, **kwargs):
25+
phys_cell = ref_to_phys.phys_cell
26+
ref_cell = ref_to_phys.ref_cell
2727
phys_element = element(phys_cell, *args, **kwargs).fiat_equivalent
2828
finat_element = element(ref_cell, *args, **kwargs)
2929

@@ -65,9 +65,8 @@ def check_zany_mapping(element, ref_cell, phys_cell, *args, **kwargs):
6565
# Zany map the results
6666
num_bfs = phys_element.space_dimension()
6767
num_dofs = finat_element.space_dimension()
68-
mappng = MyMapping(ref_cell, phys_cell)
6968
try:
70-
Mgem = finat_element.basis_transformation(mappng)
69+
Mgem = finat_element.basis_transformation(ref_to_phys)
7170
M = evaluate([Mgem])[0].arr
7271
ref_vals_zany = np.tensordot(M, ref_vals_piola, (-1, 0))
7372
except AttributeError:
@@ -86,54 +85,37 @@ def check_zany_mapping(element, ref_cell, phys_cell, *args, **kwargs):
8685
assert np.allclose(ref_vals_zany, phys_vals[:num_dofs])
8786

8887

89-
@pytest.fixture
90-
def ref_el(request):
91-
K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
92-
return K
93-
94-
95-
@pytest.fixture
96-
def phys_el(request):
97-
K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
98-
K[2].vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84))
99-
K[3].vertices = ((0, 0, 0),
100-
(1., 0.1, -0.37),
101-
(0.01, 0.987, -.23),
102-
(-0.1, -0.2, 1.38))
103-
return K
104-
105-
10688
@pytest.mark.parametrize("element", [
10789
finat.Morley,
10890
finat.Hermite,
10991
finat.Bell,
11092
])
111-
def test_C1_elements(ref_el, phys_el, element):
112-
check_zany_mapping(element, ref_el[2], phys_el[2])
93+
def test_C1_elements(ref_to_phys, element):
94+
check_zany_mapping(element, ref_to_phys[2])
11395

11496

11597
@pytest.mark.parametrize("element", [
11698
finat.QuadraticPowellSabin6,
11799
finat.QuadraticPowellSabin12,
118100
finat.ReducedHsiehCloughTocher,
119101
])
120-
def test_C1_macroelements(ref_el, phys_el, element):
102+
def test_C1_macroelements(ref_to_phys, element):
121103
kwargs = {}
122104
if element == finat.QuadraticPowellSabin12:
123105
kwargs = dict(avg=True)
124-
check_zany_mapping(element, ref_el[2], phys_el[2], **kwargs)
106+
check_zany_mapping(element, ref_to_phys[2], **kwargs)
125107

126108

127109
@pytest.mark.parametrize("element, degree", [
128110
*((finat.Argyris, k) for k in range(5, 8)),
129111
*((finat.HsiehCloughTocher, k) for k in range(3, 6))
130112
])
131-
def test_high_order_C1_elements(ref_el, phys_el, element, degree):
132-
check_zany_mapping(element, ref_el[2], phys_el[2], degree, avg=True)
113+
def test_high_order_C1_elements(ref_to_phys, element, degree):
114+
check_zany_mapping(element, ref_to_phys[2], degree, avg=True)
133115

134116

135-
def test_argyris_point(ref_el, phys_el):
136-
check_zany_mapping(finat.Argyris, ref_el[2], phys_el[2], variant="point")
117+
def test_argyris_point(ref_to_phys):
118+
check_zany_mapping(finat.Argyris, ref_to_phys[2], variant="point")
137119

138120

139121
zany_piola_elements = {
@@ -162,15 +144,15 @@ def test_argyris_point(ref_el, phys_el):
162144
*((2, e) for e in zany_piola_elements[3]),
163145
*((3, e) for e in zany_piola_elements[3]),
164146
])
165-
def test_piola(ref_el, phys_el, element, dimension):
166-
check_zany_mapping(element, ref_el[dimension], phys_el[dimension])
147+
def test_piola(ref_to_phys, element, dimension):
148+
check_zany_mapping(element, ref_to_phys[dimension])
167149

168150

169151
@pytest.mark.parametrize("element, degree, variant", [
170152
*((finat.HuZhang, k, v) for v in ("integral", "point") for k in range(3, 6)),
171153
])
172-
def test_piola_triangle_high_order(ref_el, phys_el, element, degree, variant):
173-
check_zany_mapping(element, ref_el[2], phys_el[2], degree, variant)
154+
def test_piola_triangle_high_order(ref_to_phys, element, degree, variant):
155+
check_zany_mapping(element, ref_to_phys[2], degree, variant)
174156

175157

176158
@pytest.mark.parametrize("element, degree", [
@@ -180,5 +162,5 @@ def test_piola_triangle_high_order(ref_el, phys_el, element, degree, variant):
180162
*((finat.GopalakrishnanLedererSchoberlSecondKind, k) for k in range(0, 3)),
181163
])
182164
@pytest.mark.parametrize("dimension", [2, 3])
183-
def test_affine(ref_el, phys_el, element, degree, dimension):
184-
check_zany_mapping(element, ref_el[dimension], phys_el[dimension], degree)
165+
def test_affine(ref_to_phys, element, degree, dimension):
166+
check_zany_mapping(element, ref_to_phys[dimension], degree)

0 commit comments

Comments
 (0)