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