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
36
37 from sfc.symbolic_utils import grad, ddx, symbols, symbolic_matrix
38 from sfc.geometry import UFCCell
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
44 return reduce(operator.__mul__, sequence, 1)
45
46
47
48
49
51 print type(polygon)
52 print dir(polygon)
53 print polygon.no_space_dim()
54 print polygon.str()
55
57 sfc_debug("Entering create_syfi_polygon")
58 if cell == "interval": p = SyFi.ReferenceLine()
59 elif cell == "triangle": p = SyFi.ReferenceTriangle()
60 elif cell == "tetrahedron": p = SyFi.ReferenceTetrahedron()
61 elif cell == "quadrilateral": p = SyFi.ReferenceRectangle()
62 elif cell == "hexahedron": p = SyFi.ReferenceBox()
63 else: raise ValueError("Unknown element cell '%s'." % cell)
64 sfc_debug("Leaving create_syfi_polygon")
65 return p
66
68 "Create a basic element with SyFi."
69 sfc_debug("Entering create_syfi_element")
70
71 sfc_assert(not isinstance(e, ufl.MixedElement), "Only creating SyFi elements for basic elements.")
72 f = e.family()
73
74
75 d = e.degree()
76 if not isinstance(d, int): d = default_order
77
78 if f in ("Lagrange", "CG"):
79 fe = SyFi.Lagrange(polygon, d)
80 elif f in ("Discontinuous Lagrange", "DG"):
81 if d == 0: fe = SyFi.P0(polygon, 0)
82 else: fe = SyFi.DiscontinuousLagrange(polygon, d)
83 elif f in ("Crouzeix-Raviart", "CR"):
84 fe = SyFi.CrouzeixRaviart(polygon, d)
85 elif f in ("Bubble", "B"):
86 fe = SyFi.Bubble(polygon, d)
87 elif f in ("Brezzi-Douglas-Marini", "BDM"):
88 raise NotImplementedError("Not implemented element family '%s'." % f)
89 elif f in ("Brezzi-Douglas-Fortin-Marini", "BDFM"):
90 raise NotImplementedError("Not implemented element family '%s'." % f)
91 elif f in ("Raviart-Thomas", "RT"):
92 raise NotImplementedError("Not implemented element family '%s'." % f)
93 elif f in ("Nedelec 1st kind H(div)", "N1div"):
94 raise NotImplementedError("Not implemented element family '%s'." % f)
95 elif f in ("Nedelec 2nd kind H(div)", "N2div"):
96 raise NotImplementedError("Not implemented element family '%s'." % f)
97 elif f in ("Nedelec 1st kind H(curl)", "N1curl"):
98 raise NotImplementedError("Not implemented element family '%s'." % f)
99 elif f in ("Nedelec 2nd kind H(curl)", "N2curl"):
100 raise NotImplementedError("Not implemented element family '%s'." % f)
101 elif f in ("Quadrature", "Q"):
102 raise NotImplementedError("Not implemented element family '%s'." % f)
103 elif f in ("Boundary Quadrature", "BQ"):
104 raise NotImplementedError("Not implemented element family '%s'." % f)
105 else:
106 raise NotImplementedError("Unknown element family '%s'." % f)
107
108 sfc_debug("Leaving create_syfi_element")
109 return fe
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
129 __slots__ = (
130 "options", "signature",
131 "dof_map_classname", "finite_element_classname",
132
133 "ufl_element", "syfi_element",
134
135 "quad_rule", "ufl_cell", "polygon", "cell",
136
137 "local_dimension", "geometric_dimension", "topological_dimension",
138
139 "value_shape", "value_rank", "value_size", "value_components",
140
141 "sub_elements", "sub_element_dof_offsets", "sub_element_value_offsets",
142
143 "_basis_function_cache", "_basis_function_derivative_cache",
144
145 "dof_xi", "dof_x",
146
147 "entity_dofs", "dof_entity", "num_entity_dofs", "facet_dofs", "num_facet_dofs",
148
149 "p0", "p", "G", "GinvT",
150 )
151
152 - def __init__(self, ufl_element, quad_rule=None, options=None, cache=None):
153 sfc_debug("Entering ElementRepresentation.__init__")
154
155
156 assert isinstance(ufl_element, ufl.FiniteElementBase)
157 self.ufl_element = ufl_element
158 self.quad_rule = quad_rule
159 if options is None:
160 self.options = default_parameters()
161 else:
162 self.options = options
163 if cache is None:
164 cache = {}
165
166
167 self.signature = repr(self.ufl_element)
168 self.dof_map_classname = dof_map_classname(self.ufl_element)
169 self.finite_element_classname = finite_element_classname(self.ufl_element)
170
171
172 self.ufl_cell = self.ufl_element.cell()
173 self.polygon = create_syfi_polygon(self.ufl_cell.domain())
174 self.cell = UFCCell(self.polygon)
175
176 self.geometric_dimension = self.cell.nsd
177 self.topological_dimension = self.cell.nsd
178
179
180 self.value_shape = self.ufl_element.value_shape()
181 self.value_rank = len(self.value_shape)
182 self.value_size = product(self.value_shape)
183 self.value_components = ufl.permutation.compute_indices(self.value_shape)
184
185
186 self.sub_elements = []
187 self.sub_element_dof_offsets = []
188 self.sub_element_value_offsets = []
189
190 if isinstance(self.ufl_element, ufl.MixedElement):
191 dof_offset = 0
192 value_size_offset = 0
193
194
195 for se in ufl_element.sub_elements():
196 rep = cache.get(se, None)
197 if rep is None:
198 rep = ElementRepresentation(se, self.quad_rule, self.options, cache)
199 cache[se] = rep
200 self.sub_elements.append(rep)
201
202
203 self.sub_element_dof_offsets.append(dof_offset)
204 self.sub_element_value_offsets.append(value_size_offset)
205
206 dof_offset += rep.local_dimension
207 value_size_offset += rep.value_size
208
209
210 self.sub_element_dof_offsets.append(dof_offset)
211 self.sub_element_value_offsets.append(value_size_offset)
212
213
214 self.syfi_element = None
215 self.local_dimension = dof_offset
216
217 elif self.ufl_element.family() == "Quadrature":
218
219 self.syfi_element = None
220 self.local_dimension = self.quad_rule.num_points
221
222 elif self.ufl_element.family() == "Boundary Quadrature":
223
224 sfc_error("Boundary Quadrature elements not implemented!")
225 self.syfi_element = None
226 self.local_dimension = self.facet_quad_rule.num_points
227
228 elif self.ufl_element.family() == "Real":
229
230 self.syfi_element = None
231 self.local_dimension = self.value_size
232
233 else:
234
235 self.syfi_element = create_syfi_element(self.ufl_element, self.polygon, default_order=self.options.code.finite_element.default_order_of_element)
236 self.local_dimension = self.syfi_element.nbf()
237
238
239 self._def_symbols()
240
241
242 self._precomp_coords()
243
244
245 self._build_entity_dofs()
246 self._build_facet_dofs()
247
248
249 self._basis_function_cache = {}
250 self._basis_function_derivative_cache = {}
251
252 sfc_debug("Leaving ElementRepresentation.__init__")
253
255 nsd = self.cell.nsd
256
257
258 self.p0 = swiginac.matrix(nsd, 1, symbols(["x0", "y0", "z0"][:nsd]))
259 self.p = swiginac.matrix(nsd, 1, symbols(["x", "y", "z"][:nsd]))
260
261
262
263
264 self.G = symbolic_matrix(nsd, nsd, "G")
265 self.GinvT = symbolic_matrix(nsd, nsd, "GinvT")
266
268 "Compute local dof coordinate of dof i."
269 nsd = self.cell.nsd
270
271 if self.sub_elements:
272 sub_element_index = self.dof_to_sub_element_index(i)
273 sub_dof = i - self.sub_element_dof_offsets[sub_element_index]
274 xi = self.sub_elements[sub_element_index]._dof_xi(sub_dof)
275
276 elif self.ufl_element.family() == "Quadrature":
277 xi = swiginac.matrix(nsd, 1, self.quad_rule.points[i])
278
279 elif self.ufl_element.family() == "Real":
280 xi = swiginac.matrix(nsd, 1, [0.0]*nsd)
281
282 else:
283
284 dofi = self.syfi_element.dof(i)
285 if isinstance(dofi[0], swiginac.numeric):
286
287 dof_xi_list = dofi
288 elif isinstance(dofi[0], list):
289
290 if isinstance(dofi[0][0], list):
291
292 midpoint = [0 for i in range(len(dofi[0][0]))]
293 for d in dofi[0][0:]:
294 for p, dp in enumerate(d):
295 midpoint[p] += dp
296 for p in range(len(d)):
297 midpoint[p] /= len(dofi[0])
298 dof_xi_list = midpoint
299 else:
300
301 dof_xi_list = dofi[0]
302 xi = swiginac.matrix(nsd, 1, dof_xi_list)
303
304 return xi
305
307 "Precompute dof coordinates."
308 self.dof_xi = []
309 self.dof_x = []
310
311 for i in range(self.local_dimension):
312
313 dof_xi = self._dof_xi(i)
314 self.dof_xi.append(dof_xi)
315
316
317 dof_x = (self.G.mul(dof_xi).add(self.p0)).evalm()
318 self.dof_x.append(dof_x)
319
321 "Build dof vs mesh entity relations."
322
323
324
325
326 self.entity_dofs = []
327 for i in range(self.topological_dimension+1):
328 lists = [[] for j in range(self.cell.num_entities[i])]
329 self.entity_dofs.append(lists)
330
331 if self.ufl_element.family() in ("Discontinuous Lagrange", "Quadrature"):
332
333 self.entity_dofs[self.topological_dimension][0] = list(range(self.local_dimension))
334
335 elif self.ufl_element.family() == "Boundary Quadrature":
336 sfc_error("Boundary Quadrature not handled.")
337
338 elif self.ufl_element.family() == "Real":
339 pass
340
341 else:
342
343
344 for k in range(self.local_dimension):
345 (d, i) = self.cell.find_entity(self.dof_xi[k])
346 self.entity_dofs[d][i].append(k)
347
348
349 self.dof_entity = [None]*self.local_dimension
350 for d in range(self.topological_dimension+1):
351 for i in range(self.cell.num_entities[d]):
352 for j, k in enumerate(self.entity_dofs[d][i]):
353 sfc_assert(self.dof_entity[k] is None, "Expected to set each dof only once.")
354 self.dof_entity[k] = (d, i, j)
355
356
357 self.num_entity_dofs = tuple(len(self.entity_dofs[d][0]) for d in range(self.topological_dimension+1))
358
359
360
361 for d in range(self.topological_dimension+1):
362 for doflist in self.entity_dofs[d]:
363
364
365 assert len(doflist) == self.num_entity_dofs[d]
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 N = self._basis_function_cache.get((i,component), None)
416 if N is not None:
417 return N
418
419 if self.sub_elements:
420
421 sub_element_index = self.dof_to_sub_element_index(i)
422 sub_element = self.sub_elements[sub_element_index]
423
424
425 sub_dof = i - self.sub_element_dof_offsets[sub_element_index]
426
427
428 comp_index = component_to_index(component, self.value_shape)
429 value_offset = self.sub_element_value_offsets[sub_element_index]
430
431
432 is_nonzero = (comp_index >= value_offset) and (comp_index < (value_offset + sub_element.value_size))
433 if is_nonzero:
434
435 sub_comp_index = comp_index - value_offset
436 sub_component = index_to_component(sub_comp_index, sub_element.value_shape)
437
438
439 N = sub_element.basis_function(sub_dof, sub_component)
440
441 else:
442 N = swiginac.numeric(0.0)
443
444 elif "Quadrature" in self.ufl_element.family():
445 sfc_error("Cannot compute basis functions for quadrature element.")
446 N = swiginac.numeric(1.0)
447
448 elif "Real" in self.ufl_element.family():
449
450 N = swiginac.numeric(1.0)
451
452 else:
453
454 N = self.syfi_element.N(i)
455 if isinstance(N, swiginac.matrix):
456 sfc_assert((int(N.rows()), int(N.cols())) == self.value_shape, "Shape mismatch")
457 return N[component]
458 sfc_assert(component == (), "Found scalar basic element, expecting no component, got %s." % repr(component))
459
460
461 self._basis_function_cache[(i,component)] = N
462 return N
463
465
466 directions = tuple(sorted(directions))
467
468
469 DN = self._basis_function_derivative_cache.get((i,component,directions), None)
470 if DN is not None:
471 return DN
472
473
474 DN = self.basis_function(i, component)
475 for j in directions:
476 DN = ddx(DN, j, self.GinvT)
477
478
479 self._basis_function_derivative_cache[(i,component,directions)] = DN
480 return DN
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504