1
2 """
3 This module contains code generation functionality for the ufc::*_integral classes.
4 """
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26 from itertools import tee
27 import swiginac
28 import ufc_utils
29
30 from sfc.common.output import sfc_debug, sfc_error, sfc_assert, sfc_warning
31 from sfc.common.utilities import indices_subset, unique
32
33 from sfc.codegeneration.codeformatting import (indent, CodeFormatter,
34 gen_symbol_declarations,
35 gen_const_token_definitions,
36 gen_token_declarations,
37 gen_token_assignments,
38 gen_token_additions,
39 gen_token_prints)
40
41 from sfc.codegeneration.support_code import (generate_sign_correction_support_code,
42 gen_A_reset_block)
43
44 from sfc.codegeneration.integralcgbase import IntegralCGBase
45
46
47 from sfc.quadrature import (gen_quadrature_rule_definition,
48 gen_quadrature_rule_definitions)
50 "Shared RegularIntegralCG code for several base classes."
53
71
73 if not self.options.enable_debug_code:
74 return ""
75 code = CodeFormatter()
76 code.begin_debug()
77 code.begin_block()
78 code += 'std::cout << std::endl;'
79 code += 'std::cout << "SFC DEBUGGING OUTPUT: " << std::endl;'
80 code += 'std::cout << "void %s::tabulate_tensor(...)" << std::endl;' % self.itgrep.classname
81 code += 'std::cout << "{" << std::endl;'
82 fr = self.itgrep.formrep
83 fd = fr.formdata
84 for k in range(fr.num_coefficients):
85 element = fd.elements[k]
86 rep = fr.element_reps[element]
87 code += 'std::cout << " w[%d][:] = [ ";' % k
88 code += "for(int j=0; j<%d; j++)" % rep.local_dimension
89 code.begin_block()
90 code += 'std::cout << w[%d][j] << ", ";' % k
91 code.end_block()
92 code += 'std::cout << " ]" << std::endl;'
93 code += 'std::cout << " now computing element tensor..." << std::endl;'
94 code.end_block()
95 code.end_debug()
96 return str(code)
97
99 if not self.options.enable_debug_code:
100 return ""
101 code = CodeFormatter()
102 code.begin_debug()
103 code.begin_block()
104 code += 'std::cout << " ... done computing element tensor." << std::endl;'
105 code += "for(int k=0; k<%d; k++)" % len(self.itgrep.indices)
106 code.begin_block()
107
108 code += 'std::cout << " A[" << k << "]" << " = " << A[k] << std::endl;'
109 code.end_block()
110 code += 'std::cout << "}" << std::endl;'
111 code.end_block()
112 code.end_debug()
113 return str(code)
114
116 "Generate code for a partition of the linearized computational graph."
117 code = gen_token_assignments( self.itgrep.iter_partition(data, deps,
118 basis_functions) )
119 if code:
120 deps = ", ".join(deps)
121 if basis_functions:
122 deps += "; using basis functions %s" % str(basis_functions)
123 code = "// Partition depending on %s\n{\n%s\n}" % (deps, indent(code))
124 return code
125
127 sfc_debug("Entering RegularIntegralCG.gen_geometry_block")
128
129 code = CodeFormatter()
130
131
132
133
134
135
136
137
138
139 geometry_tokens = self.itgrep.iter_geometry_tokens()
140
141
142 if self.options.enable_debug_code:
143 geometry_tokens, debug_geometry_tokens = tee(geometry_tokens)
144
145 code += "// Geometric quantities"
146 code += gen_const_token_definitions(geometry_tokens)
147
148 if self.options.enable_debug_code:
149 code += ""
150 code.begin_debug()
151 code += gen_token_prints(debug_geometry_tokens)
152 code.end_debug()
153 code += ""
154
155
156 if self.itgrep._on_facet:
157 fr = self.itgrep.formrep
158 code += "// Geometric quantities on each facet"
159 for facet in range(fr.cell.num_facets):
160 facet_tokens = self.itgrep.iter_facet_tokens(facet)
161
162 if facet == 0:
163
164 if self.options.enable_debug_code:
165 facet_tokens, debug_facet_tokens = tee(facet_tokens)
166
167
168 facet_tokens, facet_tokens2 = tee(facet_tokens)
169 code += gen_token_declarations(facet_tokens2)
170 code.begin_switch("facet")
171
172 code.begin_case(facet, braces=True)
173 code += "// Geometric quantities on facet %d" % facet
174 code += gen_token_assignments(facet_tokens)
175 code.end_case()
176 code += "default:"
177 code.indent()
178 code += 'throw std::runtime_error("Invalid facet number.");'
179 code.dedent()
180 code.end_switch()
181
182 if self.options.enable_debug_code:
183 code += ""
184 code.begin_debug()
185 code += 'std::cout << "facet = " << facet << std::endl;'
186 code += gen_token_prints(debug_facet_tokens)
187 code.end_debug()
188 code += ""
189
190 sfc_debug("Leaving RegularIntegralCG.gen_geometry_block")
191 return str(code)
192
198
200 "Assigning values to element tensor (for pre-integrated tensor components)"
201
202 assert data.integral is self.itgrep.symbolic_integral
203
204 if self.itgrep._on_facet:
205
206 fr = self.itgrep.formrep
207 code = CodeFormatter()
208 code.begin_switch("facet")
209 for facet in range(fr.cell.num_facets):
210 A_tokens = self.itgrep.iter_A_tokens(data, facet)
211 code.begin_case(facet, braces=True)
212 code += "// Integrated element tensor entries"
213 code += gen_token_assignments(A_tokens)
214 code.end_case()
215 code.end_switch()
216 code = str(code)
217
218 else:
219 A_tokens = self.itgrep.iter_A_tokens(data)
220 code = gen_token_assignments(A_tokens)
221 code = "// Integrated element tensor entries\n{\n%s\n}" % code
222
223 return code
224
226 "Beginning of quadrature loop."
227
228 assert data.integral in self.itgrep.quadrature_integrals
229
230 code = CodeFormatter()
231
232 if self.itgrep._on_facet:
233 num_points = data.facet_quad_rules[0].num_points
234 code += "for(int iq=0; iq<%d; iq++)" % num_points
235 code += "{"
236 code.indent()
237 code += "// Fetch quadrature rule on facet"
238 code += "const double quad_weight = quadrature_weights[facet][iq];"
239 code += "const double *qp = quadrature_points[facet][iq];"
240 else:
241 num_points = data.quad_rule.num_points
242 code += "for(int iq=0; iq<%d; iq++)" % num_points
243 code += "{"
244 code.indent()
245 code += "// Fetch quadrature rule"
246 code += "const double quad_weight = quadrature_weights[iq];"
247 code += "const double *qp = quadrature_points[iq];"
248
249 nsd = self.itgrep.formrep.cell.nsd
250 gr = self.itgrep.formrep.geomrep
251 if 0:
252 coords = ', '.join('qp[%d]' % i for i in range(nsd))
253 code += "const double xi[%d] = { %s };" % (nsd, coords)
254 assert gr.xi_sym[0].printc() == "xi[0]"
255 else:
256 coords = ', '.join('%s=qp[%d]' % ("xyz"[i], i) for i in range(nsd))
257 code += "const double %s;" % coords
258 assert gr.xi_sym[0].printc() == "x"
259
260 if 0:
261 coords = ', '.join(gr.x_expr[i] for i in range(nsd))
262 code += "const double x[%d] = { %s };" % (nsd, coords)
263 assert gr.x_sym[0].printc() == "x[0]"
264 elif 0:
265 coords = ', '.join("%s=%s" % (gr.x_sym[i].printc(), gr.x_expr[i].printc()) for i in range(nsd))
266 code += "const double %s;" % (coords,)
267 assert gr.x_sym[0].printc() == "x0"
268
269 code.dedent()
270
271 return str(code)
272
274 "Accumulation of nonzero element tensor entries."
275 assert data.integral in self.itgrep.quadrature_integrals
276 A_tokens = self.itgrep.iter_A_tokens(data)
277 code = gen_token_additions((A_sym, A_expr)
278 for (A_sym, A_expr) in A_tokens
279 if A_expr != 0)
280 assert code
281 code = "// Accumulating element tensor entries\n{\n%s\n}" % indent(code)
282 return code
283
286
288 if self.itgrep._symbol_counter > 0:
289 code = "double s[%d];" % self.itgrep._symbol_counter
290 else:
291 code = ""
292 return code
293
295 "Generic function to generate tabulate_tensor composed by reusable code blocks."
296
297
298
299
300
301
302
303 sfc_debug("Entering RegularIntegralCG.tabulate_tensor")
304
305 fr = self.itgrep.formrep
306 fd = fr.formdata
307 r = fr.rank
308
309 def generate_partition_blocks(data, bf_dep_list, known):
310 pblocks = []
311 known = frozenset(known)
312
313
314 todo = []
315 for (iota, keep) in bf_dep_list:
316 deps = known
317 if keep is not None:
318 deps |= frozenset("v%d" % j for (j,k) in enumerate(keep) if k)
319 if deps in data.partitions:
320 todo.append((deps, iota))
321
322
323 for t in todo:
324 deps, iota = t
325 code = self.gen_partition_block(data, deps, iota)
326 if code:
327 pblocks += [code]
328 return pblocks
329
330
331 blocks = []
332
333
334 blocks += [self.gen_pre_debug_code()]
335
336
337 blocks += [self.gen_geometry_block()]
338
339
340 integral = self.itgrep.symbolic_integral
341 if integral is not None:
342 data = self.itgrep.integral_data[integral.measure()]
343
344
345
346
347
348 blocks += [self.gen_A_assignment_block(data)]
349 else:
350 blocks += [gen_A_reset_block(len(self.itgrep.indices))]
351
352
353 bf_dep_list = [((), None)]
354
355 if r == 1:
356 m = fr.element_reps[fd.elements[0]].local_dimension
357 bf_dep_list += [((i,), (True,)) for i in range(m)]
358 elif r == 2:
359 m = fr.element_reps[fd.elements[0]].local_dimension
360 n = fr.element_reps[fd.elements[1]].local_dimension
361 bf_dep_list += [((None, i), (False, True)) for i in range(m)]
362 bf_dep_list += [((i, None), (True, False)) for i in range(m)]
363 bf_dep_list += [((i, j), (True, True)) for i in range(m)
364 for j in range(n)]
365 elif r > 2:
366 sfc_error("Support for higher order tensors not implemented.")
367
368
369 for integral in self.itgrep.quadrature_integrals:
370
371
372 data = self.itgrep.integral_data[integral.measure()]
373
374
375 blocks += generate_partition_blocks(data, bf_dep_list, ())
376
377
378
379 blocks += generate_partition_blocks(data, bf_dep_list, ("c",))
380
381
382 blocks += [self.gen_quadrature_begin_block(data)]
383
384
385 qblocks = []
386
387 if True:
388 qblocks += [indent(self.gen_quadrature_runtime_block(data))]
389
390
391
392 qblocks += generate_partition_blocks(data, bf_dep_list, ("x",))
393
394
395 qblocks += generate_partition_blocks(data, bf_dep_list, ("c", "x"))
396
397
398 qblocks += [self.gen_A_accumulation_block(data)]
399
400 blocks += [indent(code) for code in qblocks]
401
402
403 blocks += [self.gen_quadrature_end_block()]
404
405
406 blocks += [self.gen_post_debug_code()]
407
408
409 code = self.gen_symbol_allocation_block()
410 if code:
411 blocks.insert(0, code)
412
413
414 final_code = "\n\n".join(blocks)
415
416 sfc_debug("Leaving RegularIntegralCG.tabulate_tensor")
417 return final_code
418
420 return 'throw std::runtime_error("Not implemented.");'
421
422
425 RegularIntegralCG.__init__(self, itgrep)
426
427
428
429
430
431
432 for integral in self.itgrep.quadrature_integrals:
433 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!")
434 data = self.itgrep.integral_data[integral.measure()]
435
436 a, b = tee(self.itgrep.iter_member_quad_tokens(data))
437 self._iter_member_quad_tokens1 = a
438 self._iter_member_quad_tokens2 = b
439
441 return ufc_utils.cell_integral_header
442
444 return ufc_utils.cell_integral_implementation
445
447 sfc_debug("entering CellIntegral.gen_members")
448 code = ""
449
450
451
452 for integral in self.itgrep.quadrature_integrals:
453 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!")
454 data = self.itgrep.integral_data[integral.measure()]
455
456
457 member_quad_tokens = self._iter_member_quad_tokens1
458 member_quad_token_names = (str(s[0]).split(".")[-1] for s in member_quad_tokens)
459 decl = gen_symbol_declarations(member_quad_token_names)
460 if decl:
461 code += "\n\n"
462 code += "struct member_quad_tokens\n{\n"
463 code += indent(decl)
464 code += "\n};\n"
465 code += "member_quad_tokens qt[%d];\n" % data.quad_rule.num_points
466
467 sfc_debug("leaving CellIntegral.gen_members")
468 return code
469
471 sfc_debug("entering CellIntegral.gen_constructor")
472 code = ""
473
474
475
476 for integral in self.itgrep.quadrature_integrals:
477 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!")
478 data = self.itgrep.integral_data[integral.measure()]
479
480
481 member_quad_tokens = self._iter_member_quad_tokens2
482 member_quad_tokens_assignments = gen_token_assignments(member_quad_tokens)
483 if member_quad_tokens_assignments:
484
485 quad_loop_start_code = self.gen_quadrature_begin_block(data)
486
487
488 code += "// Quad loop start code:\n"
489 code += quad_loop_start_code + "\n"
490 code += "// Member quad tokens:\n"
491 code += indent(member_quad_tokens_assignments)
492 code += "\n}\n"
493
494 sfc_debug("leaving CellIntegral.gen_constructor")
495 return code
496
499
500
503 RegularIntegralCG.__init__(self, itgrep)
504
505
506
507
508
509 for integral in self.itgrep.quadrature_integrals:
510 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!")
511 data = self.itgrep.integral_data[integral.measure()]
512
513 a, b = tee(self.itgrep.iter_member_quad_tokens(data))
514 self._iter_member_quad_tokens1 = a
515 self._iter_member_quad_tokens2 = b
516
518 return ufc_utils.exterior_facet_integral_header
519
521 return ufc_utils.exterior_facet_integral_implementation
522
524 sfc_debug("entering ExteriorFacetIntegral.gen_members")
525 code = ""
526
527
528 for integral in self.itgrep.quadrature_integrals:
529 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!")
530 data = self.itgrep.integral_data[integral.measure()]
531
532 member_quad_tokens = self._iter_member_quad_tokens1
533 member_quad_token_names = (str(s[0]).split(".")[-1] for s in member_quad_tokens)
534 decl = gen_symbol_declarations(member_quad_token_names)
535 if decl:
536 code += "\n\n"
537 code += "struct member_quad_tokens\n{\n"
538 code += indent(decl)
539 code += "\n};\n"
540 code += "member_quad_tokens qt[%d][%d];\n" % (self.itgrep.formrep.cell.num_facets, data.facet_quad_rules[0].num_points)
541
542 sfc_debug("leaving ExteriorFacetIntegral.gen_members")
543 return code
544
546 sfc_debug("entering ExteriorFacetIntegral.gen_constructor")
547 code = ""
548
549
550 for integral in self.itgrep.quadrature_integrals:
551 sfc_assert(len(self.itgrep.quadrature_integrals) <= 1, "Not implemented!")
552 data = self.itgrep.integral_data[integral.measure()]
553
554 member_quad_tokens = self._iter_member_quad_tokens2
555 member_quad_tokens_assignments = gen_token_assignments(member_quad_tokens)
556 if member_quad_tokens_assignments:
557
558 facet_loop_begin = "for(int facet=0; facet<%d; facet++)\n{\n" % self.itgrep.formrep.cell.num_facets
559 quad_loop_start_code = self.gen_quadrature_begin_block(data)
560
561
562 code += facet_loop_begin
563 code += quad_loop_start_code
564 code += indent(member_quad_tokens_assignments)
565 code += "\n}\n"
566 code += "\n}\n"
567
568 sfc_debug("leaving ExteriorFacetIntegral.gen_constructor")
569 return code
570
573
574
578
580 return ufc_utils.exterior_facet_integral_header
581
583 return ufc_utils.exterior_facet_integral_implementation
584
587
590
593
595 return 'throw std::runtime_error("tabulate_tensor for interior facet integrals not implemented");'
596
598 return 'throw std::runtime_error("tabulate_tensor for interior facet integrals not implemented");'
599
600
601
602 CellIntegralCG = RegularCellIntegralCG
603 ExteriorFacetIntegralCG = RegularExteriorFacetIntegralCG
604 InteriorFacetIntegralCG = RegularInteriorFacetIntegralCG
605