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 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 shape = ""
70 if polygon.str().find("Triangle") >= 0:
71 shape = "triangle"
72 if polygon.str().find("Line") >= 0:
73 shape = "interval"
74 if polygon.str().find("Tetrahedron") >= 0:
75 shape = "tetrahedron"
76 table = find_quadrature_rule(shape, order, family="gauss")
77
78 x, y, z = symbols(["x", "y", "z"])
79 xi0, xi1, xi2 = symbols(["xi0", "xi1", "xi2"])
80
81 G, x0 = geometry_mapping(polygon)
82 xi = swiginac.matrix(nsd, 1, [xi0, xi1, xi2])
83 global_point = G*xi + x0
84 ss = function.subs(x == global_point[0,0])
85 if global_point.rows() > 1: ss = ss.subs(y == global_point[1,0])
86 if global_point.rows() > 2: ss = ss.subs(z == global_point[2,0])
87
88 sum = 0.0
89 for k in range(table.num_points):
90 p = table.points[k]
91 w = table.weights[k]
92 s = ss.subs(xi0 == p[0])
93 if len(p) > 1: s = s.subs(xi1 == p[1])
94 if len(p) > 2: s = s.subs(xi2 == p[2])
95 sum += s*w
96 sum *= G.determinant()
97 return sum
98
99 -def integral(polygon, function, integration_mode=True, integration_order=None):
100 """
101 A small wrapper on top of the functions
102 symbolic_integral and numeric_integral.
103 """
104
105 if integration_mode == "symbolic":
106 return symbolic_integral(polygon, function)
107 elif integration_mode == "quadrature":
108 return numeric_integral(polygon, function, integration_order)
109 sfc_error("Not implemented integration mode %s" % integration_mode)
110