Skip to content

Commit 681dea3

Browse files
pbrubeckpbrubeck
authored
Remove entity permutations for non-Lagrange elements (#156)
* Remove entity permutations for non-Lagrange elements * Refactor FDM element IntegratedLegendre * Deal with trivial restriction * integral(q) hierarchical variants * Legendre: fix scale * Reorder FDM basis by increasing eigenvalue
1 parent db127a4 commit 681dea3

File tree

5 files changed

+101
-114
lines changed

5 files changed

+101
-114
lines changed

‎FIAT/fdm_element.py‎

Lines changed: 42 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -11,10 +11,9 @@
1111

1212
from FIAT import dual_set, finite_element, functional, quadrature
1313
from FIAT.barycentric_interpolation import LagrangePolynomialSet
14-
from FIAT.hierarchical import IntegratedLegendre
15-
from FIAT.orientation_utils import make_entity_permutations_simplex
16-
from FIAT.P0 import P0Dual
14+
from FIAT.polynomial_set import ONPolynomialSet
1715
from FIAT.reference_element import LINE
16+
from FIAT.P0 import P0
1817

1918

2019
def sym_eig(A, B):
@@ -40,23 +39,24 @@ def tridiag_eig(A, B):
4039
numpy.reciprocal(Z, out=Z)
4140
numpy.multiply(numpy.sqrt(Z), V, out=V)
4241
numpy.multiply(V, a[:, None], out=V)
43-
return Z, V
42+
# Reorder by increasing eigenvalue
43+
return Z[::-1], V[:, ::-1]
4444

4545

4646
class FDMDual(dual_set.DualSet):
4747
"""The dual basis for 1D elements with FDM shape functions."""
4848
def __init__(self, ref_el, degree, bc_order=1, formdegree=0, orthogonalize=False):
4949
# Define the generalized eigenproblem on a reference element
50-
embedded_degree = degree + formdegree
51-
embedded = IntegratedLegendre(ref_el, embedded_degree)
52-
edim = embedded.space_dimension()
53-
self.embedded = embedded
50+
P = ONPolynomialSet(ref_el, degree + formdegree, variant="bubble")
51+
Pdim = len(P)
52+
# Apply even / odd reordering on edge bubbles
53+
P = P.take([*range(2), *range(2, Pdim, 2), *range(3, Pdim, 2)])
54+
self.poly_set = P
5455

55-
vertices = ref_el.get_vertices()
5656
if bc_order == 1 and formdegree == 0:
57-
rule = quadrature.GaussLobattoLegendreQuadratureLineRule(ref_el, edim+1)
57+
rule = quadrature.GaussLobattoLegendreQuadratureLineRule(ref_el, Pdim+1)
5858
else:
59-
rule = quadrature.GaussLegendreQuadratureLineRule(ref_el, edim)
59+
rule = quadrature.GaussLegendreQuadratureLineRule(ref_el, Pdim)
6060
self.rule = rule
6161

6262
solve_eig = sym_eig
@@ -65,21 +65,21 @@ def __init__(self, ref_el, degree, bc_order=1, formdegree=0, orthogonalize=False
6565

6666
# Tabulate the BC nodes
6767
if bc_order == 0:
68-
C = numpy.empty((0, edim), "d")
68+
C = numpy.empty((0, Pdim), "d")
6969
else:
70-
constraints = embedded.tabulate(bc_order-1, vertices)
70+
constraints = P.tabulate(ref_el.get_vertices(), bc_order-1)
7171
C = numpy.transpose(numpy.column_stack(list(constraints.values())))
7272
bdof = slice(None, C.shape[0])
7373
idof = slice(C.shape[0], None)
7474

75-
# Tabulate the basis that splits the DOFs into interior and bcs
76-
E = numpy.eye(edim)
75+
# Coefficients of the vertex and interior modes
76+
E = numpy.eye(Pdim)
7777
E[bdof, idof] = -C[:, idof]
7878
E[bdof, :] = numpy.linalg.solve(C[:, bdof], E[bdof, :])
7979

8080
# Assemble the constrained Galerkin matrices on the reference cell
8181
k = max(1, bc_order)
82-
phi = embedded.tabulate(k, rule.get_points())
82+
phi = P.tabulate(rule.get_points(), k)
8383
wts = rule.get_weights()
8484
E0 = numpy.dot(E.T, phi[(0, )])
8585
Ek = numpy.dot(E.T, phi[(k, )])
@@ -113,34 +113,33 @@ def __init__(self, ref_el, degree, bc_order=1, formdegree=0, orthogonalize=False
113113
numpy.multiply(S, lam, out=S)
114114
basis = numpy.dot(S.T, Ek)
115115

116-
bc_nodes = []
116+
sd = ref_el.get_spatial_dimension()
117+
top = ref_el.get_topology()
118+
entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
119+
nodes = []
117120
if formdegree == 0:
118121
if orthogonalize:
119122
idof = slice(None)
120123
elif bc_order > 0:
121-
bc_nodes += [functional.PointEvaluation(ref_el, x) for x in vertices]
122-
for alpha in range(1, bc_order):
123-
bc_nodes += [functional.PointDerivative(ref_el, x, [alpha]) for x in vertices]
124+
# Vertex dofs -- jet evaluation
125+
for v in sorted(top[0]):
126+
cur = len(nodes)
127+
x, = ref_el.make_points(0, v, 0)
128+
nodes.append(functional.PointEvaluation(ref_el, x))
129+
nodes.extend(functional.PointDerivative(ref_el, x, (alpha, ))
130+
for alpha in range(1, bc_order))
131+
entity_ids[0][v].extend(range(cur, len(nodes)))
124132

125133
elif bc_order > 0:
126-
basis[bdof, :] = numpy.sqrt(1.0E0 / ref_el.volume())
134+
basis[bdof] = numpy.sqrt(1.0E0 / ref_el.volume())
127135
idof = slice(formdegree, None)
128136

129-
nodes = bc_nodes + [functional.IntegralMoment(ref_el, rule, f) for f in basis[idof]]
137+
# Interior dofs -- moments against eigenfunctions
138+
cur = len(nodes)
139+
nodes.extend(functional.IntegralMoment(ref_el, rule, f) for f in basis[idof])
140+
entity_ids[sd][0].extend(range(cur, len(nodes)))
130141

131-
if len(bc_nodes) > 0:
132-
entity_ids = {0: {0: [0], 1: [1]},
133-
1: {0: list(range(2, degree+1))}}
134-
entity_permutations = {}
135-
entity_permutations[0] = {0: {0: [0]}, 1: {0: [0]}}
136-
entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree - 1)}
137-
else:
138-
entity_ids = {0: {0: [], 1: []},
139-
1: {0: list(range(0, degree+1))}}
140-
entity_permutations = {}
141-
entity_permutations[0] = {0: {0: []}, 1: {0: []}}
142-
entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree + 1)}
143-
super().__init__(nodes, ref_el, entity_ids, entity_permutations)
142+
super().__init__(nodes, ref_el, entity_ids)
144143

145144

146145
class FDMFiniteElement(finite_element.CiarletElement):
@@ -158,16 +157,18 @@ def _bc_order(self):
158157
def _formdegree(self):
159158
pass
160159

160+
def __new__(cls, ref_el, degree):
161+
if cls._formdegree == 1 and degree == 0:
162+
return P0(ref_el)
163+
return super().__new__(cls)
164+
161165
def __init__(self, ref_el, degree):
162166
if ref_el.shape != LINE:
163167
raise ValueError("%s is only defined in one dimension." % type(self))
164-
if degree == 0:
165-
dual = P0Dual(ref_el)
166-
else:
167-
dual = FDMDual(ref_el, degree, bc_order=self._bc_order,
168-
formdegree=self._formdegree, orthogonalize=self._orthogonalize)
168+
dual = FDMDual(ref_el, degree, bc_order=self._bc_order,
169+
formdegree=self._formdegree, orthogonalize=self._orthogonalize)
169170
if self._formdegree == 0:
170-
poly_set = dual.embedded.poly_set
171+
poly_set = dual.poly_set
171172
else:
172173
lr = quadrature.GaussLegendreQuadratureLineRule(ref_el, degree+1)
173174
poly_set = LagrangePolynomialSet(ref_el, lr.get_points())

‎FIAT/hierarchical.py‎

Lines changed: 42 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -8,17 +8,19 @@
88

99
import numpy
1010

11-
from FIAT import finite_element, dual_set, functional, P0
11+
from FIAT import finite_element, dual_set, functional
1212
from FIAT.reference_element import symmetric_simplex
13-
from FIAT.orientation_utils import make_entity_permutations_simplex
1413
from FIAT.quadrature import FacetQuadratureRule
1514
from FIAT.quadrature_schemes import create_quadrature
1615
from FIAT.polynomial_set import ONPolynomialSet, make_bubbles
17-
from FIAT.check_format_variant import parse_lagrange_variant
16+
from FIAT.check_format_variant import check_format_variant
17+
from FIAT.P0 import P0
1818

1919

2020
def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None):
2121
"""Tabulate the L2-duals of the hierarchical C0 basis."""
22+
if ref_el.get_spatial_dimension() == 0:
23+
degree = 0
2224
if interpolant_deg is None:
2325
interpolant_deg = degree
2426
Q = create_quadrature(ref_el, degree + interpolant_deg)
@@ -32,105 +34,86 @@ def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None):
3234

3335
class LegendreDual(dual_set.DualSet):
3436
"""The dual basis for Legendre elements."""
35-
def __init__(self, ref_el, degree, codim=0):
36-
nodes = []
37-
entity_ids = {}
38-
entity_permutations = {}
39-
37+
def __init__(self, ref_el, degree, codim=0, interpolant_deg=None):
38+
if interpolant_deg is None:
39+
interpolant_deg = degree
4040
sd = ref_el.get_spatial_dimension()
4141
top = ref_el.get_topology()
42-
for dim in sorted(top):
43-
npoints = degree + 1 if dim == sd - codim else 0
44-
perms = make_entity_permutations_simplex(dim, npoints)
45-
entity_permutations[dim] = {}
46-
entity_ids[dim] = {}
47-
if npoints == 0:
48-
for entity in sorted(top[dim]):
49-
entity_ids[dim][entity] = []
50-
entity_permutations[dim][entity] = perms
51-
continue
42+
entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
43+
nodes = []
5244

53-
ref_facet = ref_el.construct_subelement(dim)
54-
poly_set = ONPolynomialSet(ref_facet, degree)
55-
Q_ref = create_quadrature(ref_facet, 2 * degree)
56-
phis = poly_set.tabulate(Q_ref.get_points())[(0,) * dim]
57-
for entity in sorted(top[dim]):
58-
cur = len(nodes)
59-
Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
60-
nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis)
61-
entity_ids[dim][entity] = list(range(cur, len(nodes)))
62-
entity_permutations[dim][entity] = perms
45+
dim = sd - codim
46+
ref_facet = ref_el.construct_subelement(dim)
47+
poly_set = ONPolynomialSet(ref_facet, degree, scale="L2 piola")
48+
Q_ref = create_quadrature(ref_facet, degree + interpolant_deg)
49+
Phis = poly_set.tabulate(Q_ref.get_points())[(0,) * dim]
50+
for entity in sorted(top[dim]):
51+
cur = len(nodes)
52+
Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
53+
# phis must transform like a d-form to undo the measure transformation
54+
scale = 1 / Q_facet.jacobian_determinant()
55+
phis = scale * Phis
56+
nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis)
57+
entity_ids[dim][entity].extend(range(cur, len(nodes)))
6358

64-
super().__init__(nodes, ref_el, entity_ids, entity_permutations)
59+
super().__init__(nodes, ref_el, entity_ids)
6560

6661

6762
class Legendre(finite_element.CiarletElement):
6863
"""Simplicial discontinuous element with Legendre polynomials."""
6964
def __new__(cls, ref_el, degree, variant=None):
7065
if degree == 0:
71-
splitting, _ = parse_lagrange_variant(variant, integral=True)
72-
if splitting is None:
66+
splitting, variant, interpolant_deg = check_format_variant(variant, degree)
67+
if splitting is None and interpolant_deg == 0:
7368
# FIXME P0 on the split requires implementing SplitSimplicialComplex.symmetry_group_size()
74-
return P0.P0(ref_el)
69+
return P0(ref_el)
7570
return super().__new__(cls)
7671

7772
def __init__(self, ref_el, degree, variant=None):
78-
splitting, _ = parse_lagrange_variant(variant, integral=True)
73+
splitting, variant, interpolant_deg = check_format_variant(variant, degree)
7974
if splitting is not None:
8075
ref_el = splitting(ref_el)
8176
poly_set = ONPolynomialSet(ref_el, degree)
82-
dual = LegendreDual(ref_el, degree)
77+
dual = LegendreDual(ref_el, degree, interpolant_deg=interpolant_deg)
8378
formdegree = ref_el.get_spatial_dimension() # n-form
8479
super().__init__(poly_set, dual, degree, formdegree)
8580

8681

8782
class IntegratedLegendreDual(dual_set.DualSet):
8883
"""The dual basis for integrated Legendre elements."""
89-
def __init__(self, ref_el, degree):
84+
def __init__(self, ref_el, degree, interpolant_deg=None):
85+
if interpolant_deg is None:
86+
interpolant_deg = degree
87+
top = ref_el.get_topology()
88+
entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
9089
nodes = []
91-
entity_ids = {}
92-
entity_permutations = {}
9390

94-
top = ref_el.get_topology()
9591
for dim in sorted(top):
96-
perms = make_entity_permutations_simplex(dim, degree - dim)
97-
entity_ids[dim] = {}
98-
entity_permutations[dim] = {}
99-
if dim == 0 or degree <= dim:
100-
for entity in sorted(top[dim]):
101-
cur = len(nodes)
102-
pts = ref_el.make_points(dim, entity, degree)
103-
nodes.extend(functional.PointEvaluation(ref_el, pt) for pt in pts)
104-
entity_ids[dim][entity] = list(range(cur, len(nodes)))
105-
entity_permutations[dim][entity] = perms
92+
if degree <= dim:
10693
continue
107-
10894
ref_facet = symmetric_simplex(dim)
109-
Q_ref, phis = make_dual_bubbles(ref_facet, degree)
95+
Q_ref, Phis = make_dual_bubbles(ref_facet, degree, interpolant_deg=interpolant_deg)
11096
for entity in sorted(top[dim]):
11197
cur = len(nodes)
11298
Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
113-
11499
# phis must transform like a d-form to undo the measure transformation
115100
scale = 1 / Q_facet.jacobian_determinant()
116-
Jphis = scale * phis
117-
118-
nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in Jphis)
119-
entity_ids[dim][entity] = list(range(cur, len(nodes)))
120-
entity_permutations[dim][entity] = perms
101+
phis = scale * Phis
102+
nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis)
103+
entity_ids[dim][entity].extend(range(cur, len(nodes)))
121104

122-
super().__init__(nodes, ref_el, entity_ids, entity_permutations)
105+
super().__init__(nodes, ref_el, entity_ids)
123106

124107

125108
class IntegratedLegendre(finite_element.CiarletElement):
126109
"""Simplicial continuous element with integrated Legendre polynomials."""
127110
def __init__(self, ref_el, degree, variant=None):
128-
splitting, _ = parse_lagrange_variant(variant, integral=True)
111+
splitting, variant, interpolant_deg = check_format_variant(variant, degree)
129112
if splitting is not None:
130113
ref_el = splitting(ref_el)
131114
if degree < 1:
132115
raise ValueError(f"{type(self).__name__} elements only valid for k >= 1")
133116
poly_set = ONPolynomialSet(ref_el, degree, variant="bubble")
134-
dual = IntegratedLegendreDual(ref_el, degree)
117+
dual = IntegratedLegendreDual(ref_el, degree, interpolant_deg=interpolant_deg)
135118
formdegree = 0 # 0-form
136119
super().__init__(poly_set, dual, degree, formdegree)

‎FIAT/polynomial_set.py‎

Lines changed: 5 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -291,21 +291,15 @@ def make_bubbles(ref_el, degree, codim=0, shape=(), scale="L2 piola"):
291291
"""Construct a polynomial set with codim bubbles up to the given degree.
292292
"""
293293
poly_set = ONPolynomialSet(ref_el, degree, shape=shape, scale=scale, variant="bubble")
294+
if ref_el.get_spatial_dimension() == 0:
295+
return poly_set
296+
294297
entity_ids = expansions.polynomial_entity_ids(ref_el, degree, continuity="C0")
295298
sd = ref_el.get_spatial_dimension()
296299
dim = sd - codim
297-
if dim == 1:
298-
# Apply even / odd reordering on edge bubbles
299-
indices = []
300-
for entity in entity_ids[dim]:
301-
ids = entity_ids[dim][entity]
302-
indices.extend(ids[::2])
303-
indices.extend(ids[1::2])
304-
else:
305-
indices = list(chain(*entity_ids[dim].values()))
306-
300+
indices = list(chain(*entity_ids[dim].values()))
307301
if shape != ():
308-
ncomp = numpy.prod(shape)
302+
ncomp = numpy.prod(shape, dtype=int)
309303
dimPk = poly_set.get_num_members() // ncomp
310304
indices = list((numpy.array(indices)[:, None] + dimPk * numpy.arange(ncomp)[None, :]).flat)
311305
poly_set = poly_set.take(indices)

‎FIAT/restricted.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ def __init__(self, dual, indices):
2424
for dof in dofs if dof in indices]
2525
nodes = [nodes_old[i] for i in indices]
2626
self._dual = dual
27+
2728
super().__init__(nodes, ref_el, entity_ids)
2829

2930
def get_indices(self, restriction_domain, take_closure=True):
@@ -69,4 +70,4 @@ def __init__(self, element, indices=None, restriction_domain=None, take_closure=
6970
assert all(e_mapping == mapping_new[0] for e_mapping in mapping_new)
7071

7172
# Call constructor of CiarletElement
72-
super().__init__(poly_set, dual, 0, element.get_formdegree(), mapping_new[0])
73+
super().__init__(poly_set, dual, element.degree(), element.get_formdegree(), mapping_new[0])

‎finat/restricted.py‎

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -32,11 +32,19 @@ def restrict(element, domain, take_closure):
3232
@restrict.register(FiatElement)
3333
def restrict_fiat(element, domain, take_closure):
3434
try:
35-
return FiatElement(FIAT.RestrictedElement(element._element,
36-
restriction_domain=domain, take_closure=take_closure))
35+
re = FIAT.RestrictedElement(element._element,
36+
restriction_domain=domain,
37+
take_closure=take_closure)
3738
except ValueError:
3839
return null_element
3940

41+
if element.space_dimension() == re.space_dimension():
42+
# FIAT.RestrictedElement wipes out entity_permuations.
43+
# In case the restriction is trivial we return the original element
44+
# to avoid reconstructing the space with an undesired permutation.
45+
return element
46+
return FiatElement(re)
47+
4048

4149
@restrict.register(PhysicallyMappedElement)
4250
def restrict_physically_mapped(element, domain, take_closure):

0 commit comments

Comments
 (0)