Skip to content

Commit ac856a0

Browse files
pbrubeckpbrubeck
committed
Allow nodes to be passed as a dict
1 parent 0a91fb5 commit ac856a0

6 files changed

Lines changed: 77 additions & 115 deletions

File tree

‎FIAT/brezzi_douglas_marini.py‎

Lines changed: 9 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -13,39 +13,29 @@
1313

1414
class BDMDualSet(dual_set.DualSet):
1515
def __init__(self, ref_el, degree, variant, interpolant_deg):
16-
nodes = []
1716
sd = ref_el.get_spatial_dimension()
1817
top = ref_el.get_topology()
1918

20-
entity_ids = {}
2119
# set to empty
22-
for dim in top:
23-
entity_ids[dim] = {}
24-
for entity in top[dim]:
25-
entity_ids[dim][entity] = []
20+
nodes = {dim: {entity: [] for entity in top[dim]} for dim in top}
2621

2722
if variant == "integral":
2823
facet = ref_el.construct_subelement(sd-1)
2924
# Facet nodes are \int_F v\cdot n p ds where p \in P_{q}
3025
# degree is q
3126
Q_ref = create_quadrature(facet, interpolant_deg + degree)
3227
Pq = polynomial_set.ONPolynomialSet(facet, degree)
33-
Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
34-
ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq_at_qpts)
28+
ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq)
3529
for f in top[sd - 1]:
36-
cur = len(nodes)
37-
nodes.extend(ells.nodes(sd-1, f))
38-
entity_ids[sd - 1][f] = list(range(cur, len(nodes)))
30+
nodes[sd-1][f] = [ells]
3931

4032
elif variant == "point":
4133
# Define each functional for the dual set
4234
# codimension 1 facets
4335
for f in top[sd - 1]:
44-
cur = len(nodes)
4536
pts_cur = ref_el.make_points(sd - 1, f, sd + degree)
46-
nodes.extend(functional.PointScaledNormalEvaluation(ref_el, f, pt)
47-
for pt in pts_cur)
48-
entity_ids[sd - 1][f] = list(range(cur, len(nodes)))
37+
nodes[sd-1][f].extend(functional.PointScaledNormalEvaluation(ref_el, f, pt)
38+
for pt in pts_cur)
4939

5040
# internal nodes
5141
if degree > 1:
@@ -56,14 +46,12 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
5646
Q_ref = create_quadrature(cell, interpolant_deg + degree - 1)
5747
Nedel = nedelec.Nedelec(cell, degree - 1, variant)
5848
mapping, = set(Nedel.mapping())
59-
Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd]
60-
ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Ned_at_qpts, mapping=mapping)
49+
Phi = Nedel.get_nodal_basis()
50+
ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Phi, mapping=mapping)
6151
for entity in top[sd]:
62-
cur = len(nodes)
63-
nodes.extend(ells.nodes(sd, entity))
64-
entity_ids[sd][entity] = list(range(cur, len(nodes)))
52+
nodes[sd][entity] = [ells]
6553

66-
super().__init__(nodes, ref_el, entity_ids)
54+
super().__init__(nodes, ref_el)
6755

6856

6957
class BrezziDouglasMarini(finite_element.CiarletElement):

‎FIAT/dual_set.py‎

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,22 @@
1212

1313

1414
class DualSet(object):
15-
def __init__(self, nodes, ref_el, entity_ids, entity_permutations=None):
15+
def __init__(self, nodes, ref_el, entity_ids=None, entity_permutations=None):
16+
if isinstance(nodes, dict):
17+
assert entity_ids is None
18+
# flatten nodes
19+
top = ref_el.get_topology()
20+
entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
21+
entity_nodes = nodes
22+
nodes = []
23+
for dim in entity_nodes:
24+
for entity in entity_nodes[dim]:
25+
cur = len(nodes)
26+
ells = entity_nodes[dim][entity]
27+
for ell in ells:
28+
nodes.extend(ell.nodes(dim, entity))
29+
entity_ids[dim][entity].extend(range(cur, len(nodes)))
30+
1631
nodes, ref_el, entity_ids, entity_permutations = merge_entities(nodes, ref_el, entity_ids, entity_permutations)
1732
self.nodes = nodes
1833
self.ref_el = ref_el

‎FIAT/functional.py‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -141,9 +141,10 @@ def __init__(self, ref_el):
141141

142142

143143
class FacetIntegralMomentBlock(FunctionalBlock):
144-
def __init__(self, ref_el, Q_ref, Phis, mapping="L2 piola"):
144+
def __init__(self, ref_el, Q_ref, P, mapping="L2 piola"):
145+
dim = P.ref_el.get_spatial_dimension()
145146
self.Q_ref = Q_ref
146-
self.Phis = Phis
147+
self.Phis = P.tabulate(Q_ref.get_points())[(0,) * dim]
147148
self.mapping = mapping
148149
super().__init__(ref_el)
149150

@@ -216,6 +217,9 @@ def __init__(self, ref_el, target_shape, pt_dict, deriv_dict,
216217
else:
217218
self.max_deriv_order = 0
218219

220+
def nodes(self, entity_dim, entity_id):
221+
yield self
222+
219223
def evaluate(self, f):
220224
"""Obsolete and broken functional evaluation.
221225

‎FIAT/nedelec.py‎

Lines changed: 10 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -104,12 +104,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
104104
sd = ref_el.get_spatial_dimension()
105105
top = ref_el.get_topology()
106106

107-
entity_ids = {}
108107
# set to empty
109-
for dim in top:
110-
entity_ids[dim] = {}
111-
for entity in top[dim]:
112-
entity_ids[dim][entity] = []
108+
nodes = {dim: {entity: [] for entity in top[dim]} for dim in top}
113109

114110
if variant == "integral":
115111
# edge nodes are \int_F v\cdot t p ds where p \in P_{q-1}(edge)
@@ -124,32 +120,23 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
124120
facet = ref_el.construct_subelement(dim)
125121
Q_ref = create_quadrature(facet, interpolant_deg + phi_deg)
126122
Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,))
127-
Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * dim]
128123
mapping = "contravariant piola"
129-
ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Phis, mapping=mapping)
124+
ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Pqmd, mapping=mapping)
130125
for entity in top[dim]:
131-
cur = len(nodes)
132-
nodes.extend(ells.nodes(dim, entity))
133-
entity_ids[dim][entity] = list(range(cur, len(nodes)))
126+
nodes[dim][entity].append(ells)
134127

135128
elif variant == "point":
136129
for i in top[1]:
137-
cur = len(nodes)
138130
# points to specify P_k on each edge
139131
pts_cur = ref_el.make_points(1, i, degree + 1)
140-
nodes.extend(functional.PointEdgeTangentEvaluation(ref_el, i, pt)
141-
for pt in pts_cur)
142-
entity_ids[1][i] = list(range(cur, len(nodes)))
132+
nodes[1][i].extend(functional.PointEdgeTangentEvaluation(ref_el, i, pt)
133+
for pt in pts_cur)
143134

144135
if sd > 2 and degree > 1: # face tangents
145136
for i in top[2]: # loop over faces
146-
cur = len(nodes)
147137
pts_cur = ref_el.make_points(2, i, degree + 1)
148-
nodes.extend(functional.PointFaceTangentEvaluation(ref_el, i, k, pt)
149-
for k in range(2) # loop over tangents
150-
for pt in pts_cur # loop over points
151-
)
152-
entity_ids[2][i] = list(range(cur, len(nodes)))
138+
nodes[2][1].extend(functional.PointFaceTangentEvaluation(ref_el, i, k, pt)
139+
for k in range(2) for pt in pts_cur)
153140

154141
# internal nodes. These are \int_T v \cdot p dx where p \in P_{q-d}^3(T)
155142
phi_deg = degree - sd
@@ -164,13 +151,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
164151

165152
for entity in top[sd]:
166153
Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref)
167-
cur = len(nodes)
168-
nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,))
169-
for d in range(sd)
170-
for phi in Phis)
171-
entity_ids[sd][entity] = list(range(cur, len(nodes)))
154+
nodes[sd][entity].extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,))
155+
for d in range(sd) for phi in Phis)
172156

173-
super().__init__(nodes, ref_el, entity_ids)
157+
super().__init__(nodes, ref_el)
174158

175159

176160
class Nedelec(finite_element.CiarletElement):

‎FIAT/nedelec_second_kind.py‎

Lines changed: 26 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -47,77 +47,68 @@ class NedelecSecondKindDual(DualSet):
4747
def __init__(self, cell, degree, variant, interpolant_deg):
4848

4949
# Define degrees of freedom
50-
(dofs, ids) = self.generate_degrees_of_freedom(cell, degree, variant, interpolant_deg)
50+
nodes = self.generate_degrees_of_freedom(cell, degree, variant, interpolant_deg)
5151
# Call init of super-class
52-
super().__init__(dofs, cell, ids)
52+
super().__init__(nodes, cell)
5353

5454
def generate_degrees_of_freedom(self, cell, degree, variant, interpolant_deg):
55-
"Generate dofs and geometry-to-dof maps (ids)."
55+
"""Generate degrees of freedom for each entity."""
5656

57-
dofs = []
58-
ids = {}
57+
nodes = {}
5958

6059
# Extract spatial dimension and topology
6160
d = cell.get_spatial_dimension()
6261
assert (d in (2, 3)), "Second kind Nedelecs only implemented in 2/3D."
6362

6463
# Zero vertex-based degrees of freedom
65-
ids[0] = {i: [] for i in sorted(cell.topology[0])}
64+
nodes[0] = {i: [] for i in sorted(cell.topology[0])}
6665

6766
# (degree+1) degrees of freedom per entity of codimension 1 (edges)
68-
(edge_dofs, ids[1]) = self._generate_edge_dofs(cell, degree, 0, variant, interpolant_deg)
69-
dofs.extend(edge_dofs)
67+
nodes[1] = self._generate_edge_dofs(cell, degree, variant, interpolant_deg)
7068

7169
# Include face degrees of freedom if 3D
7270
if d == 3:
73-
face_dofs, ids[d-1] = self._generate_facet_dofs(d-1, cell, degree,
74-
len(dofs), variant, interpolant_deg)
75-
dofs.extend(face_dofs)
71+
nodes[d-1] = self._generate_facet_dofs(d-1, cell, degree,
72+
variant, interpolant_deg)
7673

7774
# Varying degrees of freedom (possibly zero) per cell
78-
cell_dofs, ids[d] = self._generate_facet_dofs(d, cell, degree, len(dofs), variant, interpolant_deg)
79-
dofs.extend(cell_dofs)
75+
nodes[d] = self._generate_facet_dofs(d, cell, degree, variant, interpolant_deg)
8076

81-
return (dofs, ids)
77+
return nodes
8278

83-
def _generate_edge_dofs(self, cell, degree, offset, variant, interpolant_deg):
79+
def _generate_edge_dofs(self, cell, degree, variant, interpolant_deg):
8480
"""Generate degrees of freedom (dofs) for entities of
8581
codimension 1 (edges)."""
8682

8783
if variant == "integral":
88-
return self._generate_facet_dofs(1, cell, degree, offset, variant, interpolant_deg)
84+
return self._generate_facet_dofs(1, cell, degree, variant, interpolant_deg)
8985

9086
# (degree+1) tangential component point evaluation degrees of
9187
# freedom per entity of codimension 1 (edges)
92-
dofs = []
93-
ids = {}
88+
top = cell.get_topology()
89+
nodes = {entity: [] for entity in top[1]}
9490
if variant == "point":
95-
for edge in range(len(cell.get_topology()[1])):
91+
for edge in top[1]:
9692

9793
# Create points for evaluation of tangential components
9894
points = cell.make_points(1, edge, degree + 2)
9995

10096
# A tangential component evaluation for each point
101-
dofs.extend(Tangent(cell, edge, point) for point in points)
97+
nodes[edge].extend(Tangent(cell, edge, point) for point in points)
10298

103-
# Associate these dofs with this edge
104-
i = len(points) * edge
105-
ids[edge] = list(range(offset + i, offset + i + len(points)))
99+
return nodes
106100

107-
return (dofs, ids)
108-
109-
def _generate_facet_dofs(self, dim, cell, degree, offset, variant, interpolant_deg):
101+
def _generate_facet_dofs(self, dim, cell, degree, variant, interpolant_deg):
110102
"""Generate degrees of freedom (dofs) for facets."""
111103

112-
# Initialize empty dofs and identifiers (ids)
113-
num_facets = len(cell.get_topology()[dim])
114-
dofs = []
115-
ids = {i: [] for i in range(num_facets)}
104+
# Initialize empty dofs
105+
top = cell.get_topology()
106+
nodes = {entity: [] for entity in top[dim]}
116107

117108
# Return empty info if not applicable
118109
rt_degree = degree - dim + 1
119110
if rt_degree < 1:
120-
return (dofs, ids)
111+
return nodes
121112

122113
if interpolant_deg is None:
123114
interpolant_deg = degree
@@ -136,21 +127,16 @@ def _generate_facet_dofs(self, dim, cell, degree, offset, variant, interpolant_d
136127
mapping, = set(RT.mapping())
137128

138129
# Evaluate basis functions at reference quadrature points
139-
Phis = Phi.tabulate(Q_ref.get_points())[(0,) * dim]
140-
ells = FacetIntegralMomentBlock(cell, Q_ref, Phis, mapping=mapping)
130+
ells = FacetIntegralMomentBlock(cell, Q_ref, Phi, mapping=mapping)
141131

142132
# Iterate over the facets
143-
for entity in range(num_facets):
144-
cur = offset + len(dofs)
133+
for entity in top[dim]:
145134
# Construct degrees of freedom as integral moments on this cell,
146135
# using the face quadrature weighted against the values
147136
# of the (physical) Raviart--Thomas'es on the face
148-
dofs.extend(ells.nodes(dim, entity))
149-
150-
# Assign identifiers (num RTs per face + previous edge dofs)
151-
ids[entity].extend(range(cur, offset + len(dofs)))
137+
nodes[entity].append(ells)
152138

153-
return (dofs, ids)
139+
return nodes
154140

155141

156142
class NedelecSecondKind(CiarletElement):

‎FIAT/raviart_thomas.py‎

Lines changed: 10 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -63,25 +63,18 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
6363
sd = ref_el.get_spatial_dimension()
6464
top = ref_el.get_topology()
6565

66-
entity_ids = {}
6766
# set to empty
68-
for dim in top:
69-
entity_ids[dim] = {}
70-
for entity in top[dim]:
71-
entity_ids[dim][entity] = []
67+
nodes = {dim: {entity: [] for entity in top[dim]} for dim in top}
7268

7369
if variant == "integral":
7470
facet = ref_el.construct_subelement(sd-1)
7571
# Facet nodes are \int_F v\cdot n p ds where p \in P_q
7672
q = degree - 1
7773
Q_ref = create_quadrature(facet, interpolant_deg + q)
7874
Pq = polynomial_set.ONPolynomialSet(facet, q if sd > 1 else 0)
79-
Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
80-
ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq_at_qpts)
75+
ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq)
8176
for f in top[sd - 1]:
82-
cur = len(nodes)
83-
nodes.extend(ells.nodes(sd-1, f))
84-
entity_ids[sd - 1][f] = list(range(cur, len(nodes)))
77+
nodes[sd - 1][f].append(ells)
8578

8679
# internal nodes. These are \int_T v \cdot p dx where p \in P_{q-1}^d
8780
if q > 0:
@@ -92,31 +85,23 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
9285

9386
for entity in top[sd]:
9487
Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref)
95-
cur = len(nodes)
96-
nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,))
97-
for d in range(sd)
98-
for phi in Pqm1_at_qpts)
99-
entity_ids[sd][entity] = list(range(cur, len(nodes)))
88+
nodes[sd][entity].extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,))
89+
for d in range(sd) for phi in Pqm1_at_qpts)
10090

10191
elif variant == "point":
10292
# codimension 1 facets
10393
for i in top[sd - 1]:
104-
cur = len(nodes)
10594
pts_cur = ref_el.make_points(sd - 1, i, sd + degree - 1)
106-
nodes.extend(functional.PointScaledNormalEvaluation(ref_el, i, pt)
107-
for pt in pts_cur)
108-
entity_ids[sd - 1][i] = list(range(cur, len(nodes)))
95+
nodes[sd - 1][i].extend(functional.PointScaledNormalEvaluation(ref_el, i, pt)
96+
for pt in pts_cur)
10997

11098
# internal nodes. Let's just use points at a lattice
11199
if degree > 1:
112-
cur = len(nodes)
113100
pts = ref_el.make_points(sd, 0, sd + degree - 1)
114-
nodes.extend(functional.ComponentPointEvaluation(ref_el, d, (sd,), pt)
115-
for d in range(sd)
116-
for pt in pts)
117-
entity_ids[sd][0] = list(range(cur, len(nodes)))
101+
nodes[sd][0].extend(functional.ComponentPointEvaluation(ref_el, d, (sd,), pt)
102+
for d in range(sd) for pt in pts)
118103

119-
super().__init__(nodes, ref_el, entity_ids)
104+
super().__init__(nodes, ref_el)
120105

121106

122107
class RaviartThomas(finite_element.CiarletElement):

0 commit comments

Comments
 (0)