1
2
3 """
4 This module contains functions and classes converting
5 from UFL to internal SyFi representations of elements.
6 """
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30 import operator
31 import ufl
32 import swiginac
33 import SyFi
34
35 from ufl.common import component_to_index, index_to_component, product
36
37 from sfc.symbolic_utils import grad, ddx, symbols, symbolic_matrix
38 from sfc.geometry import SFCCell
39 from sfc.common.names import finite_element_classname, dof_map_classname
40 from sfc.common import sfc_assert, sfc_info, sfc_warning, sfc_error, sfc_debug
41 from sfc.common.options import default_parameters
42 from sfc.representation.geometryrepresentation import create_syfi_polygon
43
45 "Create a basic element with SyFi."
46 sfc_debug("Entering create_syfi_element")
47
48 sfc_assert(not isinstance(ufl_element, ufl.MixedElement),
49 "Only creating SyFi elements for basic elements.")
50 f = ufl_element.family()
51
52
53
54 d = ufl_element.degree()
55 if not isinstance(d, int):
56 d = default_order
57
58 if f in ("Lagrange", "CG"):
59 fe = SyFi.Lagrange(syfi_polygon, d)
60 elif f in ("Discontinuous Lagrange", "DG"):
61 if d == 0:
62 fe = SyFi.P0(syfi_polygon, 0)
63 else:
64 fe = SyFi.DiscontinuousLagrange(syfi_polygon, d)
65 elif f in ("Crouzeix-Raviart", "CR"):
66 fe = SyFi.CrouzeixRaviart(syfi_polygon, d)
67 elif f in ("Bubble", "B"):
68 fe = SyFi.Bubble(syfi_polygon, d)
69
70 elif f in ("Brezzi-Douglas-Marini", "BDM"):
71 raise NotImplementedError("Not implemented element family '%s'." % f)
72 elif f in ("Brezzi-Douglas-Fortin-Marini", "BDFM"):
73 raise NotImplementedError("Not implemented element family '%s'." % f)
74 elif f in ("Raviart-Thomas", "RT"):
75 raise NotImplementedError("Not implemented element family '%s'." % f)
76 elif f in ("Nedelec 1st kind H(div)", "N1div"):
77 raise NotImplementedError("Not implemented element family '%s'." % f)
78 elif f in ("Nedelec 2nd kind H(div)", "N2div"):
79 raise NotImplementedError("Not implemented element family '%s'." % f)
80 elif f in ("Nedelec 1st kind H(curl)", "N1curl"):
81 raise NotImplementedError("Not implemented element family '%s'." % f)
82 elif f in ("Nedelec 2nd kind H(curl)", "N2curl"):
83 raise NotImplementedError("Not implemented element family '%s'." % f)
84 elif f in ("Quadrature", "Q"):
85 raise NotImplementedError("Not implemented element family '%s'." % f)
86 elif f in ("Boundary Quadrature", "BQ"):
87 raise NotImplementedError("Not implemented element family '%s'." % f)
88
89 else:
90 raise NotImplementedError("Unknown element family '%s'." % f)
91
92 sfc_debug("Leaving create_syfi_element")
93 return fe
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
113 __slots__ = (
114 "options", "signature",
115 "dof_map_classname", "finite_element_classname",
116
117 "ufl_element", "syfi_element",
118
119 "geomrep",
120 "quad_rule", "ufl_cell", "polygon", "cell",
121
122 "local_dimension", "geometric_dimension", "topological_dimension",
123
124 "value_shape", "value_rank", "value_size", "value_components",
125
126 "sub_elements", "sub_element_dof_offsets", "sub_element_value_offsets",
127
128 "_basis_function_cache", "_basis_function_derivative_cache",
129
130 "dof_xi", "dof_x",
131
132 "entity_dofs", "dof_entity", "num_entity_dofs", "facet_dofs", "num_facet_dofs",
133
134 "p0", "p", "G", "GinvT",
135 )
136
137 - def __init__(self, ufl_element, geomrep, quad_rule=None, options=None, cache=None):
138 sfc_debug("Entering ElementRepresentation.__init__")
139
140
141 assert isinstance(ufl_element, ufl.FiniteElementBase)
142 self.ufl_element = ufl_element
143 self.quad_rule = quad_rule
144 if options is None:
145 self.options = default_parameters()
146 else:
147 self.options = options
148 if cache is None:
149 cache = {}
150
151
152 self.signature = repr(self.ufl_element)
153 self.dof_map_classname = dof_map_classname(self.ufl_element)
154 self.finite_element_classname = finite_element_classname(self.ufl_element)
155
156
157
158 self.geomrep = geomrep
159
160 self.polygon = geomrep.syfi_polygon
161 self.cell = geomrep.sfc_cell
162
163 self.geometric_dimension = self.cell.nsd
164 self.topological_dimension = self.cell.nsd
165
166
167 self.value_shape = self.ufl_element.value_shape()
168 self.value_rank = len(self.value_shape)
169 self.value_size = product(self.value_shape)
170 self.value_components = ufl.permutation.compute_indices(self.value_shape)
171
172
173 self.sub_elements = []
174 self.sub_element_dof_offsets = []
175 self.sub_element_value_offsets = []
176
177
178 if isinstance(self.ufl_element, ufl.MixedElement):
179 dof_offset = 0
180 value_size_offset = 0
181
182
183 for se in ufl_element.sub_elements():
184 rep = cache.get(se, None)
185 if rep is None:
186 rep = ElementRepresentation(se, geomrep, self.quad_rule, self.options, cache)
187 cache[se] = rep
188 self.sub_elements.append(rep)
189
190
191 self.sub_element_dof_offsets.append(dof_offset)
192 self.sub_element_value_offsets.append(value_size_offset)
193
194 dof_offset += rep.local_dimension
195 value_size_offset += rep.value_size
196
197
198 self.sub_element_dof_offsets.append(dof_offset)
199 self.sub_element_value_offsets.append(value_size_offset)
200
201
202 self.syfi_element = None
203 self.local_dimension = dof_offset
204
205 elif self.ufl_element.family() == "Quadrature":
206
207 self.syfi_element = None
208 self.local_dimension = self.quad_rule.num_points
209
210 elif self.ufl_element.family() == "Boundary Quadrature":
211
212
213 sfc_error("Boundary Quadrature elements not implemented!")
214 self.syfi_element = None
215 self.local_dimension = self.facet_quad_rule.num_points
216
217 elif self.ufl_element.family() == "Real":
218
219 self.syfi_element = None
220 self.local_dimension = self.value_size
221
222 else:
223
224 do = self.options.code.finite_element.default_order_of_element
225 self.syfi_element = create_syfi_element(self.ufl_element, self.polygon, default_order=do)
226 self.local_dimension = self.syfi_element.nbf()
227
228
229 self._def_symbols()
230
231
232 self._precomp_coords()
233
234
235 self._build_entity_dofs()
236 self._count_entity_dofs()
237 self._build_dof_entity()
238 self._build_facet_dofs()
239
240
241 self._basis_function_cache = {}
242 self._basis_function_derivative_cache = {}
243
244 sfc_debug("Leaving ElementRepresentation.__init__")
245
247 nsd = self.cell.nsd
248
249
250
251 self.p = self.geomrep.xi_sym
252 self.p0 = swiginac.matrix(nsd, 1, symbols(["x0", "y0", "z0"][:nsd]))
253
254
255
256 self.G = self.geomrep.G_sym
257 self.GinvT = symbolic_matrix(nsd, nsd, "GinvT")
258
259
261 "Compute local dof coordinate of dof i."
262 nsd = self.cell.nsd
263
264 if self.sub_elements:
265 sub_element_index = self.dof_to_sub_element_index(i)
266 sub_dof = i - self.sub_element_dof_offsets[sub_element_index]
267 xi = self.sub_elements[sub_element_index]._dof_xi(sub_dof)
268
269 elif self.ufl_element.family() == "Quadrature":
270 xi = swiginac.matrix(nsd, 1, self.quad_rule.points[i])
271
272 elif self.ufl_element.family() == "Real":
273 xi = swiginac.matrix(nsd, 1, [0.0]*nsd)
274
275 else:
276
277 dofi = self.syfi_element.dof(i)
278 if isinstance(dofi[0], swiginac.numeric):
279
280 dof_xi_list = dofi
281 elif isinstance(dofi[0], list):
282
283 if isinstance(dofi[0][0], list):
284
285 midpoint = [0 for i in range(len(dofi[0][0]))]
286 for d in dofi[0][0:]:
287 for p, dp in enumerate(d):
288 midpoint[p] += dp
289 for p in range(len(d)):
290 midpoint[p] /= len(dofi[0])
291 dof_xi_list = midpoint
292 else:
293
294 dof_xi_list = dofi[0]
295 xi = swiginac.matrix(nsd, 1, dof_xi_list)
296
297 return xi
298
300 "Precompute dof coordinates."
301 self.dof_xi = []
302 self.dof_x = []
303
304 for i in range(self.local_dimension):
305
306 dof_xi = self._dof_xi(i)
307 self.dof_xi.append(dof_xi)
308
309
310 dof_x = (self.G.mul(dof_xi).add(self.p0)).evalm()
311 self.dof_x.append(dof_x)
312
314 "Build dof vs mesh entity relations."
315
316
317
318
319
320 entity_dofs = [ [ [] for j in range(self.cell.num_entities[i]) ]
321 for i in range(self.topological_dimension+1) ]
322
323 fam = self.ufl_element.family()
324 if fam in ("Discontinuous Lagrange", "Quadrature"):
325
326 entity_dofs[self.topological_dimension][0].extend(range(self.local_dimension))
327
328 elif fam == "Boundary Quadrature":
329 sfc_error("Boundary Quadrature not handled.")
330
331 elif fam == "Real":
332 pass
333
334 else:
335
336
337 for k in range(self.local_dimension):
338 (d, i) = self.cell.find_entity(self.dof_xi[k])
339 entity_dofs[d][i].append(k)
340
341 self.entity_dofs = entity_dofs
342
344
345 self.num_entity_dofs = tuple(len(self.entity_dofs[d][0])
346 for d in range(self.topological_dimension+1))
347
348
349
350 for d in range(self.topological_dimension+1):
351
352
353 sfc_assert(all(len(doflist) == self.num_entity_dofs[d]
354 for doflist in self.entity_dofs[d]),
355 "Entities have different number of associated dofs.")
356
358
359 dof_entity = [None]*self.local_dimension
360 for d in range(self.topological_dimension+1):
361 for i in range(self.cell.num_entities[d]):
362 for j, k in enumerate(self.entity_dofs[d][i]):
363 sfc_assert(dof_entity[k] is None, "Expected to set each dof only once.")
364 dof_entity[k] = (d, i, j)
365 self.dof_entity = dof_entity
366
368 "Build facet vs dof relations."
369
370 self.facet_dofs = [[] for j in range(self.cell.num_facets)]
371
372 if self.ufl_element.family() in ("Discontinuous Lagrange", "Quadrature", "Real"):
373 pass
374
375 elif self.ufl_element.family() == "Boundary Quadrature":
376 sfc_error("Boundary Quadrature not handled.")
377
378 else:
379
380
381 for j in range(self.cell.num_facets):
382 for (i,p) in enumerate(self.dof_xi):
383 if self.cell.facet_check(j, p):
384 self.facet_dofs[j].append(i)
385
386
387 self.num_facet_dofs = len(self.facet_dofs[0])
388
389
390 sfc_assert(all(len(fdofs) == self.num_facet_dofs for fdofs in self.facet_dofs),
391 "Not the same number of dofs on each facet. This breaks an assumption in UFC.")
392
393
394
396 "Return the index of the sub element the given dof is part of."
397 n = len(self.sub_elements)
398 sfc_assert(n, "Only mixed elements have sub elements.")
399 for i in range(n+1):
400 if dof < self.sub_element_dof_offsets[i]:
401 return i-1
402 sfc_error("Invalid dof value!")
403
405 "Return a list of all dof indices for sub element with index i."
406 sfc_assert(self.sub_elements, "Only mixed elements have sub elements.")
407 a = self.sub_element_dof_offsets[i]
408 b = self.sub_element_dof_offsets[i+1]
409 return range(a, b)
410
411
412
414
415 key = (i, component)
416 N = self._basis_function_cache.get(key)
417 if N is not None:
418 return N
419
420 if self.sub_elements:
421
422 sub_element_index = self.dof_to_sub_element_index(i)
423 sub_element = self.sub_elements[sub_element_index]
424
425
426 sub_dof = i - self.sub_element_dof_offsets[sub_element_index]
427
428
429 comp_index = component_to_index(component, self.value_shape)
430 value_offset = self.sub_element_value_offsets[sub_element_index]
431
432
433 is_nonzero = ( (comp_index >= value_offset) and
434 (comp_index < (value_offset + sub_element.value_size)) )
435 if is_nonzero:
436
437 sub_comp_index = comp_index - value_offset
438 sub_component = index_to_component(sub_comp_index, sub_element.value_shape)
439
440
441 N = sub_element.basis_function(sub_dof, sub_component)
442
443 else:
444 N = swiginac.numeric(0.0)
445
446 elif "Quadrature" in self.ufl_element.family():
447 sfc_error("Cannot compute basis functions for quadrature element.")
448 N = swiginac.numeric(1.0)
449
450 elif "Real" in self.ufl_element.family():
451
452
453 N = swiginac.numeric(1.0)
454
455 else:
456
457 N = self.syfi_element.N(i)
458 if isinstance(N, swiginac.matrix):
459 sfc_assert((int(N.rows()), int(N.cols())) == self.value_shape, "Shape mismatch")
460 return N[component]
461 sfc_assert(component == (),
462 ("Found scalar basic element, expecting no component, got %s."
463 % repr(component)))
464
465
466 self._basis_function_cache[key] = N
467 return N
468
470
471 directions = tuple(sorted(directions))
472 key = (i, component, directions)
473
474
475 DN = self._basis_function_derivative_cache.get(key, None)
476 if DN is None:
477
478 DN = self.basis_function(i, component)
479 for j in directions:
480 DN = ddx(DN, j, self.GinvT)
481
482 self._basis_function_derivative_cache[key] = DN
483
484 return DN
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508