1
2 """
3 This module contains shared code generation tools 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 from sfc.common.output import sfc_debug, sfc_error, sfc_assert, sfc_warning
28 from sfc.common.utilities import indices_subset, unique
29 from sfc.quadrature import gen_quadrature_rule_definition, gen_quadrature_rule_definitions
30 from sfc.codegeneration.codeformatting import indent, CodeFormatter, \
31 gen_const_token_definitions, gen_token_prints, \
32 gen_token_declarations, gen_token_assignments, gen_token_additions
33
38
42
46
59
83
85 raise NotImplementedError
86
89
92
94 raise NotImplementedError
95
97 raise NotImplementedError
98
100 if not self.options.enable_debug_code:
101 return ""
102 code = CodeFormatter()
103 code.begin_debug()
104 code.begin_block()
105 code += 'std::cout << std::endl;'
106 code += 'std::cout << "SFC DEBUGGING OUTPUT: " << std::endl;'
107 code += 'std::cout << "void %s::tabulate_tensor(...)" << std::endl;' % self.itgrep.classname
108 code += 'std::cout << "{" << std::endl;'
109 fr = self.itgrep.formrep
110 fd = fr.formdata
111 for k in range(fr.num_coefficients):
112 element = fd.elements[k]
113 rep = fr.element_reps[element]
114 code += 'std::cout << " w[%d][:] = [ ";' % k
115 code += "for(int j=0; j<%d; j++)" % rep.local_dimension
116 code.begin_block()
117 code += 'std::cout << w[%d][j] << ", ";' % k
118 code.end_block()
119 code += 'std::cout << " ]" << std::endl;'
120 code += 'std::cout << " now computing element tensor..." << std::endl;'
121 code.end_block()
122 code.end_debug()
123 return str(code)
124
126 if not self.options.enable_debug_code:
127 return ""
128 code = CodeFormatter()
129 code.begin_debug()
130 code.begin_block()
131 code += 'std::cout << " ... done computing element tensor." << std::endl;'
132 code += "for(int k=0; k<%d; k++)" % len(self.itgrep.indices)
133 code.begin_block()
134
135 code += 'std::cout << " A[" << k << "]" << " = " << A[k] << std::endl;'
136 code.end_block()
137 code += 'std::cout << "}" << std::endl;'
138 code.end_block()
139 code.end_debug()
140 return str(code)
141
143 "Generate code for a partition of the linearized computational graph."
144 code = gen_token_assignments( self.itgrep.iter_partition(data, deps, basis_functions) )
145 if code:
146 deps = ", ".join(deps)
147 if basis_functions:
148 deps += "; using basis functions %s" % str(basis_functions)
149 code = "// Partition depending on %s\n{\n%s\n}" % (deps, indent(code))
150 return code
151
153 sfc_debug("Entering IntegralCG.gen_geometry_block")
154
155 code = CodeFormatter()
156
157
158
159
160
161
162
163
164
165 geometry_tokens = self.itgrep.iter_geometry_tokens()
166
167
168 if self.options.enable_debug_code:
169 geometry_tokens, debug_geometry_tokens = tee(geometry_tokens)
170
171 code += "// Geometric quantities"
172 code += gen_const_token_definitions(geometry_tokens)
173
174 if self.options.enable_debug_code:
175 code += ""
176 code.begin_debug()
177 code += gen_token_prints(debug_geometry_tokens)
178 code.end_debug()
179 code += ""
180
181
182 if self.itgrep._on_facet:
183 fr = self.itgrep.formrep
184 code += "// Geometric quantities on each facet"
185 for facet in range(fr.cell.num_facets):
186 facet_tokens = self.itgrep.iter_facet_tokens(facet)
187
188 if facet == 0:
189
190 if self.options.enable_debug_code:
191 facet_tokens, debug_facet_tokens = tee(facet_tokens)
192
193
194 facet_tokens, facet_tokens2 = tee(facet_tokens)
195 code += gen_token_declarations(facet_tokens2)
196 code.begin_switch("facet")
197
198 code.begin_case(facet, braces=True)
199 code += "// Geometric quantities on facet %d" % facet
200 code += gen_token_assignments(facet_tokens)
201 code.end_case()
202 code += "default:"
203 code.indent()
204 code += 'throw std::runtime_error("Invalid facet number.");'
205 code.dedent()
206 code.end_switch()
207
208 if self.options.enable_debug_code:
209 code += ""
210 code.begin_debug()
211 code += 'std::cout << "facet = " << facet << std::endl;'
212 code += gen_token_prints(debug_facet_tokens)
213 code.end_debug()
214 code += ""
215
216 sfc_debug("Leaving IntegralCG.gen_geometry_block")
217 return str(code)
218
224
226 "Assigning values to element tensor (for pre-integrated tensor components)"
227
228 assert data.integral is self.itgrep.symbolic_integral
229
230 if self.itgrep._on_facet:
231
232 fr = self.itgrep.formrep
233 code = CodeFormatter()
234 code.begin_switch("facet")
235 for facet in range(fr.cell.num_facets):
236 A_tokens = self.itgrep.iter_A_tokens(data, facet)
237 code.begin_case(facet, braces=True)
238 code += "// Integrated element tensor entries"
239 code += gen_token_assignments(A_tokens)
240 code.end_case()
241 code.end_switch()
242 code = str(code)
243
244 else:
245 A_tokens = self.itgrep.iter_A_tokens(data)
246 code = gen_token_assignments(A_tokens)
247 code = "// Integrated element tensor entries\n{\n%s\n}" % code
248
249 return code
250
252
253 code = "// Reset element tensor\n"
254 code += "memset(A, 0, sizeof(double)*%d);" % len(self.itgrep.indices)
255 return code
256
258 "Beginning of quadrature loop."
259
260 assert data.integral in self.itgrep.quadrature_integrals
261
262 code = CodeFormatter()
263
264 if self.itgrep._on_facet:
265 num_points = data.facet_quad_rules[0].num_points
266 code += "for(int iq=0; iq<%d; iq++)" % num_points
267 code += "{"
268 code.indent()
269 code += "// Fetch quadrature rule"
270 code += "const double quad_weight = facet_quad_weights[facet][iq];"
271 code += "const double *_p = facet_quad_points[facet][iq];"
272 else:
273 num_points = data.quad_rule.num_points
274 code += "for(int iq=0; iq<%d; iq++)" % num_points
275 code += "{"
276 code.indent()
277 code += "// Fetch quadrature rule"
278 code += "const double quad_weight = quad_weights[iq];"
279 code += "const double *_p = quad_points[iq];"
280
281 nsd = self.itgrep.formrep.cell.nsd
282 if nsd > 0: code += "const double x = _p[0];"
283 if nsd > 1: code += "const double y = _p[1];"
284 if nsd > 2: code += "const double z = _p[2];"
285
286 code.dedent()
287
288 return str(code)
289
291 "Accumulation of nonzero element tensor entries."
292 assert data.integral in self.itgrep.quadrature_integrals
293 A_tokens = self.itgrep.iter_A_tokens(data)
294 code = gen_token_additions((A_sym, A_expr) for (A_sym, A_expr) in A_tokens if A_expr != 0)
295 assert code
296 code = "// Accumulating element tensor entries\n{\n%s\n}" % indent(code)
297 return code
298
301
303 if self.itgrep._symbol_counter > 0:
304 code = "double s[%d];" % self.itgrep._symbol_counter
305 else:
306 code = ""
307 return code
308
310 "Generic function to generate tabulate_tensor composed by reusable code blocks."
311
312
313
314
315
316
317
318 sfc_debug("Entering IntegralCG.tabulate_tensor")
319
320 fr = self.itgrep.formrep
321 fd = fr.formdata
322 r = fr.rank
323
324 def generate_partition_blocks(data, bf_dep_list, known):
325 pblocks = []
326 known = frozenset(known)
327
328
329 todo = []
330 for (iota, keep) in bf_dep_list:
331 deps = known
332 if keep is not None:
333 deps |= frozenset("v%d" % j for (j,k) in enumerate(keep) if k)
334 if deps in data.partitions:
335 todo.append((deps, iota))
336
337
338 for t in todo:
339 deps, iota = t
340 code = self.gen_partition_block(data, deps, iota)
341 if code:
342 pblocks += [code]
343 return pblocks
344
345
346 blocks = []
347
348
349 blocks += [self.gen_pre_debug_code()]
350
351
352 blocks += [self.gen_geometry_block()]
353
354
355 integral = self.itgrep.symbolic_integral
356 if integral is not None:
357 data = self.itgrep.integral_data[integral.measure()]
358
359
360
361
362
363 blocks += [self.gen_A_assignment_block(data)]
364 else:
365 blocks += [self.gen_A_reset_block()]
366
367
368 bf_dep_list = [((), None)]
369
370 if r == 1:
371 m = fr.element_reps[fd.elements[0]].local_dimension
372 bf_dep_list += [((i,), (True,)) for i in range(m)]
373 elif r == 2:
374 m = fr.element_reps[fd.elements[0]].local_dimension
375 n = fr.element_reps[fd.elements[1]].local_dimension
376 bf_dep_list += [((None, i), (False, True)) for i in range(m)]
377 bf_dep_list += [((i, None), (True, False)) for i in range(m)]
378 bf_dep_list += [((i, j), (True, True)) for i in range(m) for j in range(n)]
379 elif r > 2:
380 sfc_error("Support for higher order tensors not implemented.")
381
382
383 for integral in self.itgrep.quadrature_integrals:
384
385 data = self.itgrep.integral_data[integral.measure()]
386
387
388 blocks += generate_partition_blocks(data, bf_dep_list, ())
389
390
391 blocks += generate_partition_blocks(data, bf_dep_list, ("c",))
392
393
394 blocks += [self.gen_quadrature_begin_block(data)]
395
396
397 qblocks = []
398
399 if True:
400 qblocks += [indent(self.gen_quadrature_runtime_block(data))]
401
402
403 qblocks += generate_partition_blocks(data, bf_dep_list, ("x",))
404
405
406 qblocks += generate_partition_blocks(data, bf_dep_list, ("c", "x"))
407
408
409 qblocks += [self.gen_A_accumulation_block(data)]
410
411 blocks += [indent(code) for code in qblocks]
412
413
414 blocks += [self.gen_quadrature_end_block()]
415
416
417 blocks += [self.gen_post_debug_code()]
418
419
420 code = self.gen_symbol_allocation_block()
421 if code:
422 blocks.insert(0, code)
423
424
425 final_code = "\n\n".join(blocks)
426
427 sfc_debug("Leaving IntegralCG.tabulate_tensor")
428 return final_code
429
431 return 'throw std::runtime_error("Not implemented.");'
432