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 import swiginac
27 from sfc.common import sfc_assert, sfc_debug
28 from sfc.representation.integralrepresentation import IntegralRepresentation
29
33
35 fr = self.formrep
36 nsd = fr.cell.nsd
37
38
39
40
41
42 s = fr.facet_D_sym
43 e = fr.facet_D_expr[facet]
44 yield (s, e)
45
46
47 for i in range(nsd):
48 yield (fr.n_sym[i], fr.n_expr[facet][i])
49
50
51 if self.symbolic_integral is not None:
52 yield (fr.D_sym, fr.facet_D_sym)
53
55 "Compute expression for A[iota]. Overload in subclasses!"
56
57 integrand = data.integral.integrand()
58
59 if data.integration_method == "quadrature":
60 sfc_assert(facet is None, "Not expecting facet number.")
61
62 if self.options.safemode:
63 data.evaluator.update(iota)
64 integrand = data.evaluator.visit(integrand)
65 else:
66 n = len(data.G.V())
67 try:
68 integrand = data.vertex_data_set[iota][n-1]
69 except:
70 print data.vertex_data_set
71 raise RuntimeError
72
73 D = self.formrep.D_sym
74 A = integrand * D
75
76 if self.formrep.options.output.enable_debug_prints:
77 sfc_debug("In compute_A(", iota, "):")
78 sfc_debug(" data.integral.integrand() = ", data.integral.integrand())
79 sfc_debug(" integrand = ", integrand)
80 sfc_debug(" A = ", A)
81
82 elif data.integration_method == "symbolic":
83 sfc_assert(isinstance(facet, int), "Expecting facet number.")
84
85 data.evaluator.update(iota)
86 integrand = data.evaluator.visit(integrand)
87
88 D = self.formrep.D_sym
89 polygon = self.formrep.cell.facet_polygons[facet]
90 if isinstance(polygon, swiginac.matrix):
91
92 repmap = swiginac.exmap()
93 repmap[self.formrep.xi_sym[0]] = polygon[0]
94 A = integrand.subs(repmap) * D
95 else:
96 A = polygon.integrate(integrand) * D
97
98 if self.formrep.options.output.enable_debug_prints:
99 sfc_debug("In compute_A(", iota, "):")
100 sfc_debug(" data.integral.integrand() = ", data.integral.integrand())
101 sfc_debug(" integrand = ", integrand)
102 sfc_debug(" A = ", A)
103
104 return A
105