1
2 """
3 This module contains representation classes for integrals.
4 """
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28 from itertools import izip
29
30
31 import swiginac
32
33
34 from ufl.permutation import compute_indices
35 from ufl.algorithms import (Graph, partition, expand_indices2, expand_indices,
36 strip_variables, tree_format)
37 from ufl.classes import Expr, Terminal, UtilityType, SpatialDerivative, Indexed
38
39 from sfc.common import sfc_assert, sfc_warning, sfc_debug, sfc_error
40 from sfc.common.utilities import indices_subset
41 from sfc.symbolic_utils import symbols, symbol
42 from sfc.codegeneration.codeformatting import indent, row_major_index_string
43 from sfc.representation.swiginac_eval import SwiginacEvaluator, EvaluateAsSwiginac
44
45 from sfc.representation.integralrepresentationbase import IntegralRepresentationBase
46
48
49 Vcount = [len(vi) for vi in Vin]
50
51 Pset = set(P)
52 for i in P:
53 for j in Vout[i]:
54 if j in Pset:
55 Vcount[j] -= 1
56
57 return dict((i, (Vcount[i] == 0)) for i in P)
58
60 return (e.nops() == 0)
61
62
63
64
65
68 self.itgrep = itgrep
69 self.integral = integral
70
71
72 metadata = integral.measure().metadata()
73 self.integration_method = itgrep.options.integration_method
74 self.integration_order = itgrep.options.integration_order
75
76 self.quad_rule = None
77 self.facet_quad_rules = None
78
79 self.G = None
80 self.Vdeps = None
81 self.Vin_count = None
82 self.partitions = {}
83 self.evaluator = None
84 self.evaluate = None
85
86 self.V_symbols = {}
87 self.vertex_data_set = {}
88
89 - def store(self, value, i, basis_functions):
90
91 vd = self.vertex_data_set.get(basis_functions)
92 if vd is None:
93 vd = {}
94 self.vertex_data_set[basis_functions] = vd
95 vd[i] = value
96
97 - def get_vd(self, j, basis_functions):
98 jkeep = [(("v%d" % k) in self.Vdeps[j])
99 for k in range(self.itgrep.formrep.rank)]
100 iota, = indices_subset([basis_functions], jkeep)
101 if all(i is None for i in iota):
102 iota = ()
103 vd = self.vertex_data_set[iota]
104 return vd
105
107 vd = self.get_vd(j, basis_functions)
108 return vd[j]
109
111 "Delete stored data for j."
112 vd = self.get_vd(j, basis_functions)
113
114
116 - def __init__(self, integrals, formrep, on_facet):
117 IntegralRepresentationBase.__init__(self, integrals, formrep, on_facet)
118 sfc_debug("Entering IntegralRepresentation.__init__")
119
120
121 A_shape = []
122 for i in range(self.formrep.rank):
123 element = self.formrep.formdata.elements[i]
124 rep = self.formrep.element_reps[element]
125 A_shape.append(rep.local_dimension)
126 self.A_shape = tuple(A_shape)
127
128
129
130 self.indices = compute_indices(self.A_shape)
131
132
133 Asym = symbols("A[%s]" % row_major_index_string(i, self.A_shape)
134 for i in self.indices)
135 self.A_sym = dict(izip(self.indices, Asym))
136
137
138
139 if self.options.integration_method == "symbolic":
140
141 assert len(self.integrals) == 1
142 self.symbolic_integral = self.integrals[0]
143 self.quadrature_integrals = []
144 elif self.options.integration_method == "quadrature":
145
146
147 self.symbolic_integral = None
148 self.quadrature_integrals = self.integrals
149
150
151 self._free_symbols = []
152 self._symbol_counter = 0
153 self._build_integral_data()
154
155 sfc_debug("Leaving IntegralRepresentation.__init__")
156
170
176
178 data = IntegralData(self, integral)
179
180
181 data.quad_rule = self.formrep.quad_rule
182 data.facet_quad_rules = self.formrep.facet_quad_rules
183
184 if self.options.safemode:
185
186 data.evaluator = SwiginacEvaluator(self.formrep, use_symbols=True,
187 on_facet=self._on_facet)
188
189 else:
190
191 data.evaluate = EvaluateAsSwiginac(self.formrep, self, data,
192 on_facet=self._on_facet)
193
194
195 integrand = integral.integrand()
196 integrand = strip_variables(integrand)
197 if self.options.use_expand_indices2:
198 integrand = expand_indices2(integrand)
199 else:
200 integrand = expand_indices(integrand)
201
202
203 data.G = Graph(integrand)
204
205
206 V = data.G.V()
207 n = len(V)
208 data.Vin_count = [0]*n
209 for i, vs in enumerate(data.G.Vin()):
210 for j in vs:
211 vj = V[j]
212
213 data.Vin_count[i] += 1
214
215
216
217
218
219
220
221
222
223 data.partitions, data.Vdeps = partition(data.G)
224 return data
225
226
227
229 s = ""
230 s += "IntegralRepresentation:\n"
231 s += " classname: %s\n" % self.classname
232 s += " A_shape: %s\n" % self.A_shape
233 s += " indices: %s\n" % self.indices
234 s += " A_sym: %s\n" % self.A_sym
235
236
237
238
239
240
241 s += " UFL Integrals:\n"
242 for integral in self.integrals:
243 s += indent(str(integral)) + "\n\n"
244 return s.strip("\n")
245
246
247
249 """Allocate symbol(s) for vertex j. Either gets symbol from the free
250 symbol set or creates a new symbol and increases the symbol counter."""
251 if self._free_symbols:
252 return self._free_symbols.pop()
253 s = symbol("s[%d]" % self._symbol_counter)
254 self._symbol_counter += 1
255 data.V_symbols[j] = s
256 return s
257
259 "Delete stored data for j and make its symbols available again."
260 s = data.V_symbols.get(j)
261 if s is None:
262 sfc_debug("Trying to deallocate symbols that are not allocated!")
263 return
264 self._free_symbols.append(s)
265 del data.V_symbols[j]
266
267
268
270 sfc_debug("Entering IntegralRepresentation.iter_partition")
271
272 deps = frozenset(deps)
273
274 P = data.partitions.get(deps)
275 if not P:
276 sfc_debug("Leaving IntegralRepresentation.iter_partition, empty")
277 return
278
279 data.evaluate.current_basis_function = basis_functions
280
281
282 G = data.G
283 V = G.V()
284 E = G.E()
285 Vin = G.Vin()
286 Vout = G.Vout()
287
288
289 is_local = find_locals(Vin, Vout, P)
290
291
292 for i in P:
293 v = V[i]
294
295 if isinstance(v, UtilityType):
296
297 continue
298
299 if v.shape():
300
301 if not isinstance(v, (Terminal, SpatialDerivative)):
302 print "="*30
303 print "type:", type(v)
304 print "str:", str(v)
305 print "child types:", [type(V[j]) for j in Vout[i]]
306 print "child str:"
307 print "\n".join( " vertex %d: %s" % (j, str(V[j])) for j in Vout[i] )
308 print "number of parents:", len(Vin[i])
309 if len(Vin[i]) < 5:
310 print "parent types:", [type(V[j]) for j in Vin[i]]
311
312
313 sfc_error("Expecting all indexing to have been propagated to terminals?")
314 continue
315
316 if Vin[i] and all(isinstance(V[j], SpatialDerivative) for j in Vin[i]):
317
318
319
320 if not isinstance(v, (Terminal, SpatialDerivative, Indexed)):
321 print "="*30
322 print type(v)
323 print str(v)
324 sfc_error("Expecting all indexing to have been propagated to terminals?")
325 continue
326
327 if isinstance(v, (Indexed, SpatialDerivative)):
328 ops = v.operands()
329
330
331 if not all(isinstance(o, (Expr, swiginac.basic)) for o in ops):
332 print ";"*80
333 print tree_format(v)
334 print str(v)
335 print type(ops)
336 print str(ops)
337 print repr(ops)
338 print "types:"
339 print "\n".join(str(type(o)) for o in ops)
340 print ";"*80
341
342 e = data.evaluate(v, *ops)
343
344 else:
345
346
347
348 ops = []
349 for j in Vout[i]:
350 try:
351
352 e = data.fetch_storage(j, basis_functions)
353 except:
354 print "Failed to fetch expression for vertex %d," % j
355 print " V[%d] = %s" % (j, repr(V[j]))
356 print " parent V[%d] = %s" % (i, repr(V[i]))
357 raise RuntimeError
358 ops.append(e)
359 ops = tuple(ops)
360 e = data.evaluate(v, *ops)
361
362
363
364
365
366
367 for j in Vout[i]:
368 data.Vin_count[j] -= 1
369 if False:
370
371 self.free_symbols(data, j)
372 data.free_storage(j, basis_functions)
373
374
375 if is_simple(e):
376
377
378 data.store(e, i, basis_functions)
379 else:
380 if is_local[i]:
381 pass
382
383
384 s = self.allocate(data, i)
385 data.store(s, i, basis_functions)
386
387 yield (s, e)
388
389 sfc_debug("Leaving IntegralRepresentation.iter_partition")
390
391
392
393
394
395
396
397
398
399
400
402 """Return an iterator over member tokens dependent of spatial variables.
403
404 Overload in subclasses!"""
405
406
407 assert data.integration_method == "quadrature"
408
409
410
411
412
413
414
415 fr = self.formrep
416 fd = fr.formdata
417 generated = set()
418 for iarg in range(fr.rank + fr.num_coefficients):
419 elm = fd.elements[iarg]
420 rep = fr.element_reps[elm]
421 for i in range(rep.local_dimension):
422 for component in rep.value_components:
423
424 s = fr.v_sym(iarg, i, component, self._on_facet)
425 if not (s == 0 or s in generated):
426 e = fr.v_expr(iarg, i, component)
427 t = (s, e)
428 yield t
429 generated.add(s)
430
431 for d in range(fr.cell.nsd):
432 der = (d,)
433 s = fr.dv_sym(iarg, i, component, der, self._on_facet)
434 if not (s == 0 or s in generated):
435 e = fr.dv_expr(iarg, i, component, der)
436 t = (s, e)
437 yield t
438 generated.add(s)
439
441 """Return an iterator over geometry tokens independent of spatial variables."""
442
443 fr = self.formrep
444
445
446
447
448 for (ss,ee) in zip(fr.vx_sym, fr.vx_expr):
449 for i in range(ss.nops()):
450 yield (ss.op(i), ee.op(i))
451
452
453
454
455
456
457
458 (ss,ee) = (fr.G_sym, fr.G_expr)
459 for i in range(ss.nops()):
460 yield (ss.op(i), ee.op(i))
461
462
463 yield (fr.detGtmp_sym, fr.detGtmp_expr)
464 yield (fr.detG_sym, fr.detG_expr)
465
466
467 (ss,ee) = (fr.Ginv_sym, fr.Ginv_expr)
468 for i in range(ss.nops()):
469 yield (ss.op(i), ee.op(i))
470
471 if self._on_facet:
472
473 yield (fr.detG_sign_sym, fr.detG_sign_expr)
474 else:
475 if self.symbolic_integral is not None:
476
477 yield (fr.D_sym, fr.detG_sym)
478
480 "Return an iterator over runtime tokens dependent of spatial variables. Overload in subclasses!"
481 assert data.integration_method == "quadrature"
482
483 fr = self.formrep
484 gr = fr.geomrep
485 fd = fr.formdata
486 nsd = gr.sfc_cell.nsd
487
488
489
490
491
492
493 (ss,ee) = (gr.x_sym, gr.x_expr)
494 for i in range(ss.nops()):
495 yield (ss.op(i), ee.op(i))
496
497
498 generated = set()
499 for iarg in range(fr.rank + fr.num_coefficients):
500 elm = fd.elements[iarg]
501 rep = fr.element_reps[elm]
502 for i in range(rep.local_dimension):
503 for component in rep.value_components:
504
505 for d in range(fr.cell.nsd):
506 der = (d,)
507 s = fr.Dv_sym(iarg, i, component, der, self._on_facet)
508 if not (s == 0 or s in generated):
509 e = fr.Dv_expr(iarg, i, component, der, True, self._on_facet)
510 t = (s, e)
511 yield t
512 generated.add(s)
513
514
515 generated = set()
516 for iarg in range(fr.num_coefficients):
517 elm = fd.elements[fr.rank+iarg]
518 rep = fr.element_reps[elm]
519 for component in rep.value_components:
520
521 s = fr.w_sym(iarg, component)
522 if not s in generated:
523 e = fr.w_expr(iarg, component, True, self._on_facet)
524 t = (s, e)
525 yield t
526 generated.add(s)
527
528 for d in range(fr.cell.nsd):
529 der = (d,)
530 s = fr.Dw_sym(iarg, component, der)
531 if not (s == 0 or s in generated):
532 e = fr.Dw_expr(iarg, component, der, True, self._on_facet)
533 t = (s, e)
534 yield t
535 generated.add(s)
536
537
538 if self._on_facet:
539
540
541
542 D_expr = fr.quad_weight_sym*fr.facet_D_sym
543 else:
544
545 D_expr = fr.quad_weight_sym*fr.detG_sym
546 yield (fr.D_sym, D_expr)
547
548
549
551 "Iterate over all A[iota] tokens."
552 for iota in self.indices:
553 A_sym = self.A_sym[iota]
554 A_expr = self.compute_A(data, iota, facet)
555 yield (A_sym, A_expr)
556
557 - def compute_A(self, data, iota, facet=None):
558 "Compute expression for A[iota]. Overload in subclasses!"
559 raise NotImplementedError
560
561
563 - def __init__(self, integrals, formrep):
565
566 - def compute_A(self, data, iota, facet=None):
567 "Compute expression for A[iota]. Overload in subclasses!"
568 if self.options.safemode:
569 integrand = data.integral.integrand()
570 data.evaluator.update(iota)
571 integrand = data.evaluator.visit(integrand)
572 else:
573 n = len(data.G.V())
574 integrand = data.vertex_data_set[iota][n-1]
575
576 D = self.formrep.D_sym
577 A = integrand * D
578
579 if self.formrep.options.output.enable_debug_prints:
580 sfc_debug("In compute_A(", iota, "):")
581 sfc_debug(" data.integral.integrand() = ", data.integral.integrand())
582 sfc_debug(" integrand = ", integrand)
583 sfc_debug(" A = ", A)
584
585 return A
586
588 - def __init__(self, integrals, formrep):
590
591 - def compute_A(self, data, iota, facet=None):
592 "Compute expression for A[iota]. Overload in subclasses!"
593
594 integrand = data.integral.integrand()
595 data.evaluator.update(iota)
596 integrand = data.evaluator.visit(integrand)
597 detG = self.formrep.detG_sym
598 polygon = self.formrep.cell.polygon
599 A = polygon.integrate(integrand) * detG
600
601 if self.formrep.options.output.enable_debug_prints:
602 sfc_debug("In compute_A(", iota, "):")
603 sfc_debug(" data.integral.integrand() = ", data.integral.integrand())
604 sfc_debug(" integrand = ", integrand)
605 sfc_debug(" A = ", A)
606
607 return A
608
609
611 - def __init__(self, integrals, formrep):
613
615 fr = self.formrep
616 gr = fr.geomrep
617 nsd = fr.cell.nsd
618
619
620
621
622
623 s = gr.facet_D_sym
624 e = gr.facet_D_expr[facet]
625 yield (s, e)
626
627
628 for i in range(nsd):
629 yield (gr.n_tmp_sym[i], gr.n_tmp_expr[facet][i])
630 yield (gr.nnorm_sym, gr.nnorm_expr)
631 for i in range(nsd):
632 yield (gr.n_sym[i], gr.n_expr[i])
633
634
635 if self.symbolic_integral is not None:
636 yield (gr.D_sym, gr.facet_D_sym)
637
638
640 - def __init__(self, integrals, formrep):
642
643 - def compute_A(self, data, iota, facet=None):
644 "Compute expression for A[iota]. Overload in subclasses!"
645
646 integrand = data.integral.integrand()
647
648 sfc_assert(facet is None, "Not expecting facet number.")
649
650 if self.options.safemode:
651 data.evaluator.update(iota)
652 integrand = data.evaluator.visit(integrand)
653 else:
654 n = len(data.G.V())
655 try:
656 integrand = data.vertex_data_set[iota][n-1]
657 except:
658 print data.vertex_data_set
659 raise RuntimeError
660
661 D = self.formrep.D_sym
662 A = integrand * D
663
664 if self.formrep.options.output.enable_debug_prints:
665 sfc_debug("In compute_A(", iota, "):")
666 sfc_debug(" data.integral.integrand() = ", data.integral.integrand())
667 sfc_debug(" integrand = ", integrand)
668 sfc_debug(" A = ", A)
669
670 return A
671
672
674 - def __init__(self, integrals, formrep):
676
677 - def compute_A(self, data, iota, facet=None):
678 "Compute expression for A[iota]. Overload in subclasses!"
679
680 integrand = data.integral.integrand()
681
682 sfc_assert(isinstance(facet, int), "Expecting facet number.")
683
684 data.evaluator.update(iota)
685 integrand = data.evaluator.visit(integrand)
686
687 D = self.formrep.D_sym
688 polygon = self.formrep.cell.facet_polygons[facet]
689 if isinstance(polygon, swiginac.matrix):
690
691 repmap = swiginac.exmap()
692 repmap[self.formrep.xi_sym[0]] = polygon[0]
693 A = integrand.subs(repmap) * D
694 else:
695 A = polygon.integrate(integrand) * D
696
697 if self.formrep.options.output.enable_debug_prints:
698 sfc_debug("In compute_A(", iota, "):")
699 sfc_debug(" data.integral.integrand() = ", data.integral.integrand())
700 sfc_debug(" integrand = ", integrand)
701 sfc_debug(" A = ", A)
702
703 return A
704
705
707 - def __init__(self, integrals, formrep):
709
710
718
726