1
2
3 """
4 This module contains code generation tools for the ufc::finite_element class.
5 """
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27 from itertools import izip
28 import ufl
29 import swiginac
30
31 from sfc.codegeneration.codeformatting import indent, CodeFormatter, gen_switch, gen_token_assignments, gen_const_token_definitions
32 from sfc.symbolic_utils import symbol, symbols, TempSymbolContext
33 from sfc.geometry import gen_geometry_code
34 from sfc.common import sfc_warning, sfc_assert, sfc_error, sfc_info
35 from sfc.common.utilities import unique
36 import ufc_utils
37
40
42 """Code generator for ufc::finite_element implementations."""
44 self.rep = elementrep
45 self.classname = elementrep.finite_element_classname
46 self.signature = repr(self.rep.ufl_element)
47 self.options = self.rep.options.code.finite_element
48 self.DN_order = self.options.evaluate_basis_derivatives_order
49
53
57
61
63 return ufc_utils.finite_element_header
64
66 return ufc_utils.finite_element_implementation
67
69
70 vars = {
71 "classname" : self.classname,
72 "constructor" : indent(self.gen_constructor()),
73 "constructor_arguments" : indent(self.gen_constructor_arguments()),
74 "initializer_list" : indent(self.gen_initializer_list()),
75 "destructor" : indent(self.gen_destructor()),
76 "create" : indent(self.gen_create()),
77 "signature" : indent(self.gen_signature()),
78 "cell_shape" : indent(self.gen_cell_shape()),
79 "geometric_dimension" : indent(self.gen_geometric_dimension()),
80 "topological_dimension" : indent(self.gen_topological_dimension()),
81 "space_dimension" : indent(self.gen_space_dimension()),
82 "value_rank" : indent(self.gen_value_rank()),
83 "value_dimension" : indent(self.gen_value_dimension()),
84 "map_from_reference_cell" : indent(self.gen_map_from_reference_cell()),
85 "map_to_reference_cell" : indent(self.gen_map_to_reference_cell()),
86 "evaluate_basis" : indent(self.gen_evaluate_basis()),
87 "evaluate_basis_all" : indent(self.gen_evaluate_basis_all()),
88 "evaluate_basis_derivatives" : indent(self.gen_evaluate_basis_derivatives()),
89 "evaluate_basis_derivatives_all" : indent(self.gen_evaluate_basis_derivatives_all()),
90 "evaluate_dof" : indent(self.gen_evaluate_dof()),
91 "evaluate_dofs" : indent(self.gen_evaluate_dofs()),
92 "interpolate_vertex_values" : indent(self.gen_interpolate_vertex_values()),
93 "num_sub_elements" : indent(self.gen_num_sub_elements()),
94 "create_sub_element" : indent(self.gen_create_sub_element()),
95 "members" : indent(self.gen_members()),
96 }
97 return vars
98
101
104
107
110
113
115 code = "return new %s();" % self.classname
116 return code
117
119 """const char* signature() const"""
120 return 'return "%s";' % self.signature
121
123 """shape cell_shape() const"""
124 return "return ufc::%s;" % self.rep.cell.shape
125
127 return "return %d;" % self.rep.cell.nsd
128
130 return "return %d;" % self.rep.cell.nsd
131
133 """unsigned int space_dimension() const"""
134 return "return %d;" % self.rep.local_dimension
135
137 """unsigned int value_rank() const"""
138 return "return %d;" % self.rep.value_rank
139
141 """unsigned int value_dimension(unsigned int i) const"""
142 if self.rep.value_rank == 0:
143 code = 'throw std::runtime_error("Rank 0 value has no dimension.");\n'
144
145
146 else:
147 dims = self.rep.value_shape
148 cases = [(i, "return %d;" % d) for (i,d) in enumerate(dims)]
149 code = gen_switch("i", cases)
150 code += 'throw std::runtime_error("Invalid dimension for rank %d value.");\n' % self.rep.value_rank
151
152 return code
153
155 return 'throw std::runtime_error("Not implemented.");'
156
158 return 'throw std::runtime_error("Not implemented.");'
159
161 """void evaluate_basis(unsigned int i,
162 double* values,
163 const double* coordinates,
164 const cell& c) const
165 """
166 if not self.options.enable_evaluate_basis:
167 return 'throw std::runtime_error("evaluate_basis not implemented.");'
168
169 nsd = self.rep.cell.nsd
170 nbf = self.rep.local_dimension
171 value_shape = self.rep.value_shape
172 value_size = self.rep.value_size
173
174 if self.rep.ufl_element.family() == "Real":
175 code = []
176 if value_size > 1:
177 code += ["memset(values, 0, sizeof(double)*%d);" % value_size]
178 code += ['values[i] = 1.0;']
179 return '\n'.join(code)
180
181
182 val_sym = symbols(["values[%d]" % d for d in range(value_size)])
183
184
185 code = CodeFormatter()
186
187 code += gen_geometry_code(nsd, detG=False, GinvT=True)
188
189 coordinates = swiginac.matrix(nsd, 1, symbols(["coordinates[%d]" % i for i in range(nsd)]))
190 xi = self.rep.GinvT.transpose().mul(coordinates - self.rep.p0).evalm()
191
192 for i in range(nsd):
193 code += "const double %s = %s;" % (("x","y","z")[i], xi[i].printc())
194
195 if value_size > 1:
196 code += "memset(values, 0, sizeof(double)*%d);" % value_size
197 code += ""
198
199
200 if nbf>1:
201 code.begin_switch("i")
202
203
204 for i in range(nbf):
205
206
207 values = []
208 for c in self.rep.value_components:
209 values.append( self.rep.basis_function(i, c) )
210 values_tokens = []
211 for d in range(value_size):
212 if not values[d].expand().is_zero():
213 values_tokens.append( (val_sym[d], values[d]) )
214
215
216 if self.options.optimize_basis:
217
218 temp_symbol = TempSymbolContext()
219 values_tokens1 = []
220 values_tokens2 = []
221 for s, e in values_tokens:
222 ts = temp_symbol()
223 values_tokens1.append( (ts, e) )
224 values_tokens2.append( (s, ts) )
225 values_tokens = None
226
227
228 values_tokens1, repmap = cse(values_tokens1, temp_symbol)
229
230
231 values_code = gen_const_token_definitions(values_tokens1) + "\n"
232 values_code += gen_token_assignments(values_tokens2)
233
234 else:
235 values_code = gen_token_assignments(values_tokens)
236
237
238
239 if nbf>1:
240 code.begin_case(i, braces=True)
241 code.new_line( indent(values_code) )
242 code.end_case()
243
244
245 else:
246 code += values_code
247
248 if nbf>1:
249 code.end_switch()
250
251 return str(code)
252
254 """/// Evaluate all basis functions at given point in cell
255 virtual void evaluate_basis_all(double* values,
256 const double* coordinates,
257 const cell& c) const
258 """
259
260 code = CodeFormatter()
261 code += "for(unsigned i = 0; i < %d; i++)" % self.rep.local_dimension
262 code.begin_block()
263 code += "evaluate_basis(i, values+i*%d, coordinates, c);" % self.rep.value_size
264 code.end_block()
265 return str(code)
266
268 """/// Evaluate order n derivatives of basis function i at given point in cell
269 void evaluate_basis_derivatives(unsigned int i,
270 unsigned int n,
271 double* values,
272 const double* coordinates,
273 const ufc::cell& c) const
274 """
275 if not self.options.enable_evaluate_basis_derivatives:
276 return 'throw std::runtime_error("evaluate_basis_derivatives not implemented.");'
277
278 nsd = self.rep.cell.nsd
279 nbf = self.rep.local_dimension
280 value_shape = self.rep.value_shape
281 value_size = self.rep.value_size
282
283 if self.rep.ufl_element.family() == "Real":
284 return '\n'.join('values[%d] = 0.0;' % d for d in range(value_size))
285
286 sfc_assert(self.DN_order <= 2, "Don't support computing higher order derivatives (yet).")
287
288 code = CodeFormatter()
289 code.begin_if("n > 2")
290 code += 'throw std::runtime_error("evaluate_basis_derivatives not implemented for the wanted derivative order.");'
291 code.end_if()
292
293
294 code += gen_geometry_code(nsd, detG=False, GinvT=True)
295
296
297 p = symbols(["x", "y", "z"][:nsd])
298 coords = symbols(["coordinates[%d]" % i for i in range(nsd)])
299 code += gen_const_token_definitions( izip(p, coords) )
300
301
302 code.begin_switch("n")
303
304 for order in range(1, self.DN_order+1):
305 code.begin_case(order, braces=True)
306
307
308 num_derivatives = nsd ** order
309 do_memset = (value_size*num_derivatives) > 6
310 if do_memset:
311 code += "memset(values, 0, sizeof(double) * %d * %d);" % (value_size, num_derivatives)
312 code += ""
313
314
315 code.begin_switch("i")
316 for ibf in range(nbf):
317 code.begin_case(ibf, braces=True)
318 temp_symbol = TempSymbolContext()
319
320 DN_tokens = []
321 all_directions = ufl.permutation.compute_permutations(self.DN_order, nsd)
322 for (j,directions) in enumerate(all_directions):
323 for (k,c) in enumerate(self.rep.value_components):
324 DN = self.rep.basis_function_derivative(ibf, c, directions)
325
326 if not (do_memset and DN.expand().is_zero()):
327 values_sym = symbol("values[%d * %d + %d]" % (j, value_size, k))
328 DN_tokens.append( (values_sym, DN) )
329
330 if self.options.optimize_basis:
331
332 DN_tokens1 = []
333 DN_tokens2 = []
334 for s, e in DN_tokens:
335 ts = temp_symbol()
336 DN_tokens1.append( (ts, e) )
337 DN_tokens2.append( (s, ts) )
338
339 DN_tokens1, repmap = cse(DN_tokens1, temp_symbol)
340
341 code += gen_const_token_definitions(DN_tokens1)
342 code += gen_token_assignments(DN_tokens2)
343 else:
344 code += gen_token_assignments(DN_tokens)
345 code.end_case()
346 code.end_switch()
347 code.end_case()
348 code += "default:"
349 code.indent()
350 code += 'throw std::runtime_error("Derivatives of this order are not supported in evaluate_basis_derivatives.");'
351 code.dedent()
352 code.end_switch()
353 return str(code)
354
356 """/// Evaluate order n derivatives of all basis functions at given point in cell
357 virtual void evaluate_basis_derivatives_all(unsigned int n,
358 double* values,
359 const double* coordinates,
360 const cell& c) const
361 """
362 code = CodeFormatter()
363 code.begin_switch("n")
364 for order in range(1, self.DN_order+1):
365 code.begin_case(order, braces=True)
366 offset = self.rep.value_size * (self.rep.cell.nsd ** order)
367 code += "for(unsigned i = 0; i < %d; i++)" % self.rep.local_dimension
368 code.begin_block()
369 code += "evaluate_basis_derivatives(n, i, values+i*%d, coordinates, c);" % offset
370 code.end_block()
371 code.end_case()
372 code += "default:"
373 code.indent()
374 code += 'throw std::runtime_error("Derivatives of this order are not implemented in evaluate_basis_derivatives_all.");'
375 code.dedent()
376 code.end_switch()
377 return str(code)
378
380 """double evaluate_dof(unsigned int i,
381 const function& f,
382 const cell& c) const
383 This implementation is general for all elements with point evaluation dofs.
384
385 TODO: Implement support for normal and tangential component dofs.
386 TODO: Implement for elements without point evaluation dofs.
387 """
388
389 nsd = self.rep.cell.nsd
390 nbf = self.rep.local_dimension
391
392
393 code = CodeFormatter()
394 code.new_text( gen_geometry_code(nsd, detG=False) )
395 code += "double v[%d];" % self.rep.value_size
396 code += "double x[%d];" % nsd
397
398
399 if nbf == 1:
400
401 for k in range(nsd):
402 code += "x[%d] = %s;" % (k, self.rep.dof_x[0][k].printc())
403 else:
404 code.begin_switch("i")
405 for i in range(nbf):
406 code.begin_case(i)
407
408 for k in range(nsd):
409 code += "x[%d] = %s;" % (k, self.rep.dof_x[i][k].printc())
410 code.end_case()
411 code.end_switch()
412
413
414 code += "f.evaluate(v, x, c);"
415
416
417
418 if nbf == 1:
419 code += "return v[0];"
420 else:
421 code += "return v[i / %d];" % (nbf // self.rep.value_size)
422
423 return str(code)
424
426 """/// Evaluate linear functionals for all dofs on the function f
427 virtual void evaluate_dofs(double* values,
428 const function& f,
429 const cell& c) const
430 """
431 code = CodeFormatter()
432 code += "for(unsigned i=0; i<%d; i++)" % self.rep.local_dimension
433 code.begin_block()
434 code += "values[i] = evaluate_dof(i, f, c);"
435 code.end_block()
436 return str(code)
437
439 """void interpolate_vertex_values(double* vertex_values,
440 const double* dof_values,
441 const cell& c) const
442 """
443
444 nbf = self.rep.local_dimension
445 nsd = self.rep.cell.nsd
446 nv = self.rep.cell.num_entities[0]
447 value_size = self.rep.value_size
448 value_shape = self.rep.value_shape
449
450
451 dof_values_sym = symbols("dof_values[%d]" % i for i in xrange(nbf))
452 vertex_values_sym = symbols("vertex_values[%d]" % i for i in xrange(nv*value_size))
453
454
455 p = self.rep.p
456
457
458 u = []
459 for component in self.rep.value_components:
460 u.append( sum(self.rep.basis_function(j, component)*dof_values_sym[j] for j in range(nbf)) )
461
462
463 vertex_values = []
464 repmap = swiginac.exmap()
465 for i in range(nv):
466
467 vx = self.rep.polygon.vertex(i)
468 for k in range(nsd):
469 repmap[p[k]] = vx[k]
470
471 for uval in u:
472 vertex_values.append(uval.subs(repmap))
473
474 code = gen_token_assignments( izip(vertex_values_sym, vertex_values) )
475 return code
476
478 """unsigned int num_sub_elements() const"""
479 return "return %d;" % len(self.rep.sub_elements)
480
482 """finite_element* create_sub_element(unsigned int i) const"""
483 if len(self.rep.sub_elements) > 1:
484 code = CodeFormatter()
485 code.begin_switch("i")
486 for i, fe in enumerate(self.rep.sub_elements):
487 code += "case %d: return new %s();" % (i, fe.finite_element_classname)
488 code.end_switch()
489 code += 'throw std::runtime_error("Invalid index in create_sub_element.");'
490 else:
491 code = "return new %s();" % self.classname
492 return str(code)
493
496
497
498
499
500
501
502
503