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.quadrature import find_quadrature_rule
31 from sfc.symbolic_utils import symbol, symbols
32
34 """
35 This function computes the geometry mapping from
36 the corresponding reference polygon to polygon.
37 It returns the tuple (G,x0) in the
38 mapping x = G xi + x0
39 """
40
41 vlen = SyFi.cvar.nsd
42
43
44 if isinstance(polygon, swiginac.matrix):
45
46 assert vlen == 1
47 assert polygon.rows() == 1
48 assert polygon.cols() == 1
49 G = swiginac.matrix(1, 1, [1])
50 x0 = polygon
51 return G, x0
52
53 Gl = []
54 v0 = polygon.vertex(0)
55 v1 = polygon.vertex(1)
56
57 if not vlen == len(v0):
58 print
59 print vlen
60 print v0
61 print v1
62 print polygon
63 print
64
65 assert vlen == len(v0)
66
67 if isinstance(polygon, SyFi.Line):
68 Gl.extend(v1[k] - v0[k] for k in range(vlen))
69
70 elif isinstance(polygon, SyFi.Triangle):
71 if vlen == 2 or vlen == 3:
72 v2 = polygon.vertex(2)
73 Gl.extend(v1[k] - v0[k] for k in range(vlen))
74 Gl.extend(v2[k] - v0[k] for k in range(vlen))
75 else:
76 raise RuntimeError("Unsupported vertex size.")
77
78 elif isinstance(polygon, SyFi.Rectangle):
79 if vlen == 2 or vlen == 3:
80 v2 = polygon.vertex(2)
81 v3 = polygon.vertex(3)
82 Gl.extend(v1[k] - v0[k] for k in range(vlen))
83 Gl.extend(v3[k] - v0[k] for k in range(vlen))
84 else:
85 raise RuntimeError("Unsupported vertex size.")
86
87 elif isinstance(polygon, SyFi.Tetrahedron):
88 if vlen != 3:
89 raise RuntimeError("Unsupported vertex size.")
90 v2 = polygon.vertex(2)
91 v3 = polygon.vertex(3)
92 Gl.extend(v1[k] - v0[k] for k in range(vlen))
93 Gl.extend(v2[k] - v0[k] for k in range(vlen))
94 Gl.extend(v3[k] - v0[k] for k in range(vlen))
95
96 elif isinstance(polygon, SyFi.Box):
97 if vlen != 3:
98 raise RuntimeError("Unsupported vertex size.")
99 v2 = polygon.vertex(2)
100 v3 = polygon.vertex(3)
101 v4 = polygon.vertex(4)
102 Gl.extend(v1[k] - v0[k] for k in range(vlen))
103 Gl.extend(v3[k] - v0[k] for k in range(vlen))
104 Gl.extend(v4[k] - v0[k] for k in range(vlen))
105
106 else:
107 raise RuntimeError("Unsupported polygon type.")
108
109 G = swiginac.matrix(vlen, vlen, Gl).transpose()
110 x0 = swiginac.matrix(vlen, 1, v0)
111
112
113
114 return G, x0
115