1
2 """
3 This module contains functions for computing
4 integrals both symbolically and numerically.
5 """
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27 import swiginac
28 import SyFi
29
30 from sfc.common import sfc_error
31 from sfc.quadrature import find_quadrature_rule
32 from sfc.geometry import UFCCell, geometry_mapping
33 from sfc.symbolic_utils.symbolic_utils import symbol, symbols
34
36 """
37 This function computes the integral of the specified function over
38 the specified polygon symbolically.
39 The function may take as arguments the global coordinates (x,y,z)
40 as well as coordinates on the reference polygon (xi0, xi1, xi2).
41 """
42 nsd = polygon.no_space_dim()
43 G, x0 = geometry_mapping(polygon)
44
45
46 x, y, z = symbols(["x", "y", "z"])
47 xi0, xi1, xi2 = symbols(["xi0", "xi1", "xi2"])
48 p = swiginac.matrix(nsd, 1, [x,y,z])
49 Ginv = G.inverse()
50
51 xi = Ginv*(p-x0)
52 ss = function.subs(xi0 == xi[0,0])
53 if xi.rows() > 1: ss = ss.subs(xi1 == xi[1,0])
54 if xi.rows() > 2: ss = ss.subs(xi2 == xi[2,0])
55
56 return polygon.integrate(ss).evalf()
57
59 """
60 This function computes the integral of the specified function over
61 the specified polygon numerically.
62 The function may take as arguments the global coordinates (x,y,z)
63 as well as coordinates on the reference polygon (xi0, xi1, xi2).
64 """
65
66
67
68 nsd = polygon.no_space_dim()
69
70
71
72
73 shape = ""
74 if polygon.str().find("Triangle") >= 0:
75 shape = "triangle"
76 if polygon.str().find("Line") >= 0:
77 shape = "interval"
78 if polygon.str().find("Tetrahedron") >= 0:
79 shape = "tetrahedron"
80 table = find_quadrature_rule(shape, order, family="gauss")
81
82 x, y, z = symbols(["x", "y", "z"])
83 xi0, xi1, xi2 = symbols(["xi0", "xi1", "xi2"])
84
85 G, x0 = geometry_mapping(polygon)
86 xi = swiginac.matrix(nsd, 1, [xi0, xi1, xi2])
87 global_point = G*xi + x0
88 ss = function.subs(x == global_point[0,0])
89 if global_point.rows() > 1: ss = ss.subs(y == global_point[1,0])
90 if global_point.rows() > 2: ss = ss.subs(z == global_point[2,0])
91
92 sum = 0.0
93 for k in range(table.num_points):
94 p = table.points[k]
95 w = table.weights[k]
96 s = ss.subs(xi0 == p[0])
97 if len(p) > 1: s = s.subs(xi1 == p[1])
98 if len(p) > 2: s = s.subs(xi2 == p[2])
99 sum += s*w
100 sum *= G.determinant()
101 return sum
102
103 -def integral(polygon, function, integration_mode=True, integration_order=None):
104 """
105 A small wrapper on top of the functions
106 symbolic_integral and numeric_integral.
107 """
108
109 if integration_mode == "symbolic":
110 return symbolic_integral(polygon, function)
111 elif integration_mode == "quadrature":
112 return numeric_integral(polygon, function, integration_order)
113 sfc_error("Not implemented integration mode %s" % integration_mode)
114