Skip to content
This repository was archived by the owner on Feb 21, 2022. It is now read-only.

Commit 490e6b9

Browse files
celdredwence-
celdred
authored and
wence-
committed
proper testing and EGL/Edge elements
1 parent 9251e55 commit 490e6b9

File tree

6 files changed

+391
-1
lines changed

6 files changed

+391
-1
lines changed

‎FIAT/__init__.py‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,8 @@
3838
from FIAT.mixed import MixedElement # noqa: F401
3939
from FIAT.restricted import RestrictedElement # noqa: F401
4040
from FIAT.quadrature_element import QuadratureElement # noqa: F401
41+
from FIAT.extended_gauss_legendre import ExtendedGaussLegendre
42+
from FIAT.edge import EdgeGaussLobattoLegendre, EdgeExtendedGaussLegendre
4143

4244
# Important functionality
4345
from FIAT.quadrature import make_quadrature # noqa: F401
@@ -65,6 +67,9 @@
6567
"Lagrange": Lagrange,
6668
"Gauss-Lobatto-Legendre": GaussLobattoLegendre,
6769
"Gauss-Legendre": GaussLegendre,
70+
"Extended-Gauss-Legendre": ExtendedGaussLegendre,
71+
"Gauss-Lobatto-Legendre Edge": EdgeGaussLobattoLegendre,
72+
"Extended-Gauss-Legendre Edge": EdgeExtendedGaussLegendre,
6873
"Morley": Morley,
6974
"Nedelec 1st kind H(curl)": Nedelec,
7075
"Nedelec 2nd kind H(curl)": NedelecSecondKind,

‎FIAT/edge.py‎

Lines changed: 176 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,176 @@
1+
# Copyright (C) 2019 Chris Eldred (Inria Grenoble Rhone-Alpes) and Werner Bauer (Inria Rennes)
2+
#
3+
# This file is part of FIAT.
4+
#
5+
# FIAT is free software: you can redistribute it and/or modify
6+
# it under the terms of the GNU Lesser General Public License as published by
7+
# the Free Software Foundation, either version 3 of the License, or
8+
# (at your option) any later version.
9+
#
10+
# FIAT is distributed in the hope that it will be useful,
11+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
# GNU Lesser General Public License for more details.
14+
#
15+
# You should have received a copy of the GNU Lesser General Public License
16+
# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
17+
#
18+
19+
import sympy
20+
import numpy as np
21+
from FIAT.finite_element import FiniteElement
22+
from FIAT.dual_set import make_entity_closure_ids
23+
from FIAT.polynomial_set import mis
24+
from FIAT import quadrature
25+
26+
x = sympy.symbols('x')
27+
28+
29+
def _lagrange_poly(i, xi):
30+
'''Returns the Langrange polynomial P_i(x) of pts xi
31+
which interpolates the values (0,0,..1,..,0) where 1 is at the ith support point
32+
:param i: non-zero location
33+
:param xi: set of N support points
34+
'''
35+
index = list(range(len(xi)))
36+
index.pop(i)
37+
return sympy.prod([(x-xi[j][0])/(xi[i][0]-xi[j][0]) for j in index])
38+
39+
40+
def _lagrange_basis(spts):
41+
symbas = []
42+
for i in range(len(spts)):
43+
symbas.append(_lagrange_poly(i, spts))
44+
return symbas
45+
46+
47+
def _create_compatible_l2_basis(cg_symbas):
48+
ncg_basis = len(cg_symbas)
49+
symbas = []
50+
for i in range(1, ncg_basis):
51+
basis = 0
52+
for j in range(i):
53+
basis = basis + sympy.diff(-cg_symbas[j])
54+
symbas.append(basis)
55+
return symbas
56+
57+
58+
class _EdgeElement(FiniteElement):
59+
60+
def __new__(cls, ref_el, degree):
61+
dim = ref_el.get_spatial_dimension()
62+
if dim == 1:
63+
self = super().__new__(cls)
64+
return self
65+
else:
66+
raise IndexError("only intervals supported for _IntervalElement")
67+
68+
def tabulate(self, order, points, entity=None):
69+
70+
if entity is None:
71+
entity = (self.ref_el.get_spatial_dimension(), 0)
72+
entity_dim, entity_id = entity
73+
dim = self.ref_el.get_spatial_dimension()
74+
75+
transform = self.ref_el.get_entity_transform(entity_dim, entity_id)
76+
points = list(map(transform, points))
77+
78+
phivals = {}
79+
for o in range(order + 1):
80+
alphas = mis(dim, o)
81+
for alpha in alphas:
82+
try:
83+
basis = self.basis[alpha]
84+
except KeyError:
85+
basis = sympy.diff(self.basis[(0,)], x)
86+
self.basis[alpha] = basis
87+
T = np.zeros((len(basis), len(points)))
88+
for i in range(len(points)):
89+
subs = {x: points[i][0]}
90+
for j, f in enumerate(basis):
91+
T[j, i] = f.evalf(subs=subs)
92+
phivals[alpha] = T
93+
94+
return phivals
95+
96+
def degree(self):
97+
return self._degree
98+
99+
def get_nodal_basis(self):
100+
raise NotImplementedError("get_nodal_basis not implemented for _IntervalElement")
101+
102+
def get_dual_set(self):
103+
raise NotImplementedError("get_dual_set is not implemented for _IntervalElement")
104+
105+
def get_coeffs(self):
106+
raise NotImplementedError("get_coeffs not implemented for _IntervalElement")
107+
108+
def entity_dofs(self):
109+
"""Return the map of topological entities to degrees of
110+
freedom for the finite element."""
111+
return self.entity_ids
112+
113+
def entity_closure_dofs(self):
114+
"""Return the map of topological entities to degrees of
115+
freedom on the closure of those entities for the finite element."""
116+
return self.entity_closure_ids
117+
118+
def value_shape(self):
119+
return ()
120+
121+
def dmats(self):
122+
raise NotImplementedError
123+
124+
def get_num_members(self, arg):
125+
raise NotImplementedError
126+
127+
def space_dimension(self):
128+
return len(self.basis[(0,)])
129+
130+
def __init__(self, ref_el, degree):
131+
132+
dim = ref_el.get_spatial_dimension()
133+
topology = ref_el.get_topology()
134+
135+
if not dim == 1:
136+
raise IndexError("only intervals supported for DMSE")
137+
138+
formdegree = 1
139+
140+
# This is general code to build empty entity_ids
141+
entity_ids = {}
142+
for dim in range(dim+1):
143+
entity_ids[dim] = {}
144+
for entity in sorted(topology[dim]):
145+
entity_ids[dim][entity] = []
146+
147+
# The only dofs for DMSE are with the interval!
148+
entity_ids[dim][0] = list(range(degree + 1))
149+
150+
# Build the basis
151+
# This is a dictionary basis[o] that gives the "basis" functions for derivative order o, where o=0 is just the basis itself
152+
# This is filled as needed for o > 0 by tabulate
153+
154+
self.entity_ids = entity_ids
155+
self.entity_closure_ids = make_entity_closure_ids(ref_el, entity_ids)
156+
self._degree = degree
157+
158+
super(_EdgeElement, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree)
159+
160+
161+
class EdgeGaussLobattoLegendre(_EdgeElement):
162+
def __init__(self, ref_el, degree):
163+
super(EdgeGaussLobattoLegendre, self).__init__(ref_el, degree)
164+
cont_pts = quadrature.GaussLobattoLegendreQuadratureLineRule(ref_el, degree+2).pts
165+
cont_basis = _lagrange_basis(cont_pts)
166+
basis = _create_compatible_l2_basis(cont_basis)
167+
self.basis = {(0,): sympy.Array(basis)}
168+
169+
170+
class EdgeExtendedGaussLegendre(_EdgeElement):
171+
def __init__(self, ref_el, degree):
172+
super(EdgeExtendedGaussLegendre, self).__init__(ref_el, degree)
173+
cont_pts = quadrature.ExtendedGaussLegendreQuadratureLineRule(ref_el, degree+2).pts
174+
cont_basis = _lagrange_basis(cont_pts)
175+
basis = _create_compatible_l2_basis(cont_basis)
176+
self.basis = {(0,): sympy.Array(basis)}

‎FIAT/extended_gauss_legendre.py‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
# Copyright (C) 2019 Chris Eldred (Inria Grenoble Rhone-Alpes) and Werner Bauer (Inria Rennes)
2+
#
3+
# This file is part of FIAT.
4+
#
5+
# FIAT is free software: you can redistribute it and/or modify
6+
# it under the terms of the GNU Lesser General Public License as published by
7+
# the Free Software Foundation, either version 3 of the License, or
8+
# (at your option) any later version.
9+
#
10+
# FIAT is distributed in the hope that it will be useful,
11+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
# GNU Lesser General Public License for more details.
14+
#
15+
# You should have received a copy of the GNU Lesser General Public License
16+
# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
17+
#
18+
# Written by Chris Eldred (Inria Grenoble Rhone-Alpes) and Werner Bauer (Inria Rennes) 2019
19+
20+
from FIAT import finite_element, polynomial_set, dual_set, functional, quadrature
21+
from FIAT.reference_element import LINE
22+
23+
24+
class ExtendedGaussLegendreDualSet(dual_set.DualSet):
25+
"""The dual basis for 1D continuous elements with nodes at the
26+
Extended-Gauss-Legendre points."""
27+
def __init__(self, ref_el, degree):
28+
entity_ids = {0: {0: [0], 1: [degree]},
29+
1: {0: list(range(1, degree))}}
30+
lr = quadrature.ExtendedGaussLegendreQuadratureLineRule(ref_el, degree+1)
31+
nodes = [functional.PointEvaluation(ref_el, x) for x in lr.pts]
32+
33+
super(ExtendedGaussLegendreDualSet, self).__init__(nodes, ref_el, entity_ids)
34+
35+
36+
class ExtendedGaussLegendre(finite_element.CiarletElement):
37+
"""1D continuous element with nodes at the Extended-Gauss-Legendre points."""
38+
def __init__(self, ref_el, degree):
39+
if ref_el.shape != LINE:
40+
raise ValueError("Extended-Gauss-Legendre elements are only defined in one dimension.")
41+
poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
42+
dual = ExtendedGaussLegendreDualSet(ref_el, degree)
43+
formdegree = 0 # 0-form
44+
super(ExtendedGaussLegendre, self).__init__(poly_set, dual, degree, formdegree)

‎FIAT/quadrature.py‎

Lines changed: 34 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -108,7 +108,7 @@ def __init__(self, ref_el, m):
108108

109109

110110
class GaussLegendreQuadratureLineRule(QuadratureRule):
111-
"""Produce the Gauss--Legendre quadrature rules on the interval using
111+
"""Produce the Gauss-Legendre quadrature rules on the interval using
112112
the implementation in numpy. This facilitates implementing
113113
discontinuous spectral elements.
114114
@@ -134,6 +134,39 @@ def __init__(self, ref_el, m):
134134
QuadratureRule.__init__(self, ref_el, xs, ws)
135135

136136

137+
class ExtendedGaussLegendreQuadratureLineRule(QuadratureRule):
138+
"""Produce the Extended-Gauss-Legendre quadrature rules on the interval using
139+
the implementation in numpy. This facilitates implementing
140+
dual continuous mimetic spectral elements.
141+
142+
The quadrature rule uses m points for a degree of precision of 2(m-2)-3 = 2m - 7.
143+
Thus, not recommend for actual use (but is useful to define the Extended-Gauss-Legendre elements)
144+
"""
145+
def __init__(self, ref_el, m):
146+
if m < 3:
147+
raise ValueError(
148+
"Extended-Gauss-Legendre quadrature invalid for fewer than 3 points")
149+
150+
xs_ref, ws_ref = numpy.polynomial.legendre.leggauss(m-2)
151+
A, b = reference_element.make_affine_mapping(((-1.,), (1.)),
152+
ref_el.get_vertices())
153+
154+
mapping = lambda x: numpy.dot(A, x) + b
155+
156+
scale = numpy.linalg.det(A)
157+
158+
ws = list([scale * w for w in ws_ref])
159+
xs = [tuple(mapping(x_ref)[0]) for x_ref in xs_ref]
160+
161+
xs.insert(0, (0.,))
162+
xs.append((1.,))
163+
# The integral of these basis functions over the interval is 0!
164+
ws.insert(0, 0.)
165+
ws.append(0.)
166+
167+
QuadratureRule.__init__(self, ref_el, tuple(xs), tuple(ws))
168+
169+
137170
class CollapsedQuadratureTriangleRule(QuadratureRule):
138171
"""Implements the collapsed quadrature rules defined in
139172
Karniadakis & Sherwin by mapping products of Gauss-Jacobi rules

‎test/unit/test_edge.py‎

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
# Copyright (C) 2016 Imperial College London and others
2+
#
3+
# This file is part of FIAT.
4+
#
5+
# FIAT is free software: you can redistribute it and/or modify
6+
# it under the terms of the GNU Lesser General Public License as published by
7+
# the Free Software Foundation, either version 3 of the License, or
8+
# (at your option) any later version.
9+
#
10+
# FIAT is distributed in the hope that it will be useful,
11+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
# GNU Lesser General Public License for more details.
14+
#
15+
# You should have received a copy of the GNU Lesser General Public License
16+
# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
17+
#
18+
# Authors:
19+
#
20+
# Christopher Eldred
21+
22+
import pytest
23+
import math
24+
import sympy
25+
import FIAT
26+
27+
x = sympy.symbols('x')
28+
29+
30+
@pytest.mark.parametrize("degree", range(0, 7))
31+
def test_gll_edge_basis_values(degree):
32+
"""Ensure that edge elements are histopolatory"""
33+
34+
s = FIAT.ufc_simplex(1)
35+
36+
fe = FIAT.EdgeGaussLobattoLegendre(s, degree)
37+
38+
basis = fe.basis[(0,)]
39+
cont_pts = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule(s, degree + 2).pts
40+
41+
for i in range(len(cont_pts)-1):
42+
for j in range(basis.shape[0]):
43+
int_sub = sympy.integrate(basis[j], (x, cont_pts[i][0], cont_pts[i+1][0]))
44+
if i == j:
45+
assert(math.isclose(int_sub, 1.))
46+
else:
47+
assert(math.isclose(int_sub, 0., abs_tol=1e-9))
48+
49+
50+
@pytest.mark.parametrize("degree", range(2, 7))
51+
def test_egl_edge_basis_values(degree):
52+
"""Ensure that edge elements are histopolatory"""
53+
54+
s = FIAT.ufc_simplex(1)
55+
56+
fe = FIAT.EdgeExtendedGaussLegendre(s, degree)
57+
58+
basis = fe.basis[(0,)]
59+
cont_pts = FIAT.quadrature.ExtendedGaussLegendreQuadratureLineRule(s, degree + 2).pts
60+
61+
for i in range(len(cont_pts)-1):
62+
for j in range(basis.shape[0]):
63+
int_sub = sympy.integrate(basis[j], (x, cont_pts[i][0], cont_pts[i+1][0]))
64+
if i == j:
65+
assert(math.isclose(int_sub, 1.))
66+
else:
67+
assert(math.isclose(int_sub, 0., abs_tol=1e-9))
68+
69+
70+
if __name__ == '__main__':
71+
import os
72+
pytest.main(os.path.abspath(__file__))

0 commit comments

Comments
 (0)