1
2
3 """
4 This module contains utility functions for symbolic computations,
5 mostly a thin layer on top of swiginac.
6 """
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28 import swiginac
29 import SyFi
30
31 from sfc.symbolic_utils.symbol_factory import symbol, symbols, symbolic_vector, symbolic_matrix
32
33
34
35 _p = symbols( ["x", "y", "z"] )
36
37
38
40 """Checks if the argument is a matrix or lst."""
41 return isinstance(A, list) \
42 or isinstance(A, tuple) \
43 or isinstance(A, swiginac.lst) \
44 or isinstance(A.evalm(), swiginac.matrix)
45
48
50 """Convert A to a swiginac.matrix from a lst or list."""
51 if isinstance(A, swiginac.lst): return swiginac.lst_to_matrix(A).evalm()
52 if isinstance(A, list): return swiginac.matrix(len(A), 1, A).evalm()
53 if isinstance(A, tuple): return swiginac.matrix(len(A), 1, A).evalm()
54 if isinstance(A.evalm(), swiginac.matrix): return A.evalm()
55 raise TypeError("ERROR: cannot convert A to a swiginac.matrix, A is of type %s" % str(type(A)))
56
58 """Returns a m x n matrix with zeros."""
59 return swiginac.matrix(m, n)
60
62 """Returns a m x n matrix with ones."""
63 A = swiginac.matrix(m, n)
64 for i in xrange(m):
65 for j in xrange(n):
66 A[i,j] = 1
67 return A
68
70 """Returns the n x n identity matrix."""
71 return swiginac.unit_matrix(n)
72
78
84
90
92 """Takes the cross product of two vectors."""
93 a = as_matrix(a)
94 b = as_matrix(b)
95 if len(a) != 3 or len(b) != 3:
96 raise TypeError("Need 3D vectors a and b in cross(a,b).")
97 im = [a[1], a[2],
98 b[1], b[2]]
99 jm = [a[2], a[0],
100 b[2], b[0]]
101 km = [a[0], a[1],
102 b[0], b[1]]
103 return swiginac.matrix(3, 1, [swiginac.matrix(2,2,m).determinant() for m in (im, jm, km)])
104
106 """Takes the inner product of A and B: multiplication for scalars, dot product for vectors and lsts, contraction for matrices."""
107 if is_indexed_type(A): A = as_matrix(A)
108 if is_indexed_type(B): B = as_matrix(B)
109 if isinstance(A, swiginac.matrix) and isinstance(B, swiginac.matrix):
110 if len(A) == len(B):
111 return sum( [A.op(i)*B.op(i) for i in xrange(len(A))] )
112 if A.cols() == 1: A = A.transpose()
113 if B.rows() == 1: B = B.transpose()
114 if A.cols() == B.rows():
115 return mul(A, B)
116 raise ValueError("Unmatching size of operands")
117 return (A*B).evalm()
118
120 return inner(a, b)
121
123 """See inner(A,B)."""
124 return inner(A,B)
125
127 """Returns the determinant of the argument."""
128 return A.evalm().determinant()
129
131 """Returns the absolute value of the argument."""
132 return swiginac.abs(x)
133
135 """Returns the transpose of the argument."""
136 return swiginac.transpose(A.evalm())
137
139 """Returns the trace of the argument."""
140 return swiginac.trace(A.evalm())
141
143 """Returns the inverse of the argument."""
144 A = A.evalm()
145 if isinstance(A, swiginac.matrix):
146 return A.evalm().inverse()
147 return 1.0 / A
148
150 if isinstance(x, swiginac.matrix):
151 r = x.rows()
152 c = x.cols()
153 if r > 1 and c > 1: return 2
154 if r*c > 1: return 1
155 return 0
156
158 if isinstance(x, swiginac.matrix):
159 r = x.rows()
160 c = x.cols()
161 if r > 1 and c > 1: return (r, c)
162 if r*c > 1: return (r*c,)
163 return (1,)
164
166 """Returns df/dx, where x is a symbol or symbolic matrix."""
167 if type(x) == swiginac.matrix:
168 f_rank = rank(f)
169 x_rank = rank(x)
170 if f_rank+x_rank > 2: raise RuntimeError("Cannot apply diff to f,x with ranks %d,%d" % (f_rank, x_rank))
171
172 A = zeros(x.rows(), x.cols())
173 for i in range(A.rows()):
174 for j in range(A.cols()):
175 sym = x[i,j]
176 assert isinstance(sym, swiginac.symbol)
177 val = f
178 A[i,j] = swiginac.diff(val, sym)
179 return A
180 assert isinstance(x, swiginac.symbol)
181 return swiginac.diff(f, x)
182
183 -def ddx(f, i, GinvT=None):
184 """Returns df/dx_i, where i is the number of the coordinate.
185 GinvT is an optional geometry mapping."""
186
187 if GinvT is None:
188 return diff(f, x(i))
189
190 dfdx = 0
191 for k in range(GinvT.rows()):
192 dfdx += GinvT[i, k] * swiginac.diff(f, _p[k])
193 return dfdx
194
195
196 -def grad(u, GinvT=None):
197 """Returns the gradient of u w.r.t. the spacial variables x,y[,z], depending on the dimension specified by SyFi.initSyFi(nsd)."""
198
199 ran = rank(u)
200 sh = shape(u)
201
202 if ran == 0:
203 n = SyFi.cvar.nsd
204 r, c = n, 1
205 u = swiginac.matrix(1,1,[u])
206 elif ran == 1:
207 n = sh[0]
208 r, c = n, n
209 else:
210 raise ValueError("Taking gradient of a matrix is not supported.")
211
212 if GinvT is None:
213 GinvT = Id(n)
214
215 g = zeros(r, c)
216 for i in range(r):
217 for j in range(c):
218 f = u[j,0]
219 for k in range(n):
220 g[i,j] += GinvT[i,k] * diff(f, _p[k])
221 return g
222
223
224 -def div(u, GinvT=None):
225 r = u.rows()
226 c = u.cols()
227 ran = rank(u)
228 if ran == 2:
229 if r != c:
230 raise RuntimeError("Taking divergence of an unsymmetric matrix.")
231 n = r
232 if GinvT is None:
233 GinvT = Id(n)
234 d = zeros(r, 1)
235 for j in range(r):
236 for i in range(r):
237 for k in range(r):
238 d[j,0] += GinvT[i,k] * diff(u[i,j], _p[k])
239 elif ran == 1:
240 n = r*c
241 if GinvT is None:
242 GinvT = Id(n)
243 d = 0
244 for i in range(n):
245 for k in range(n):
246 d += GinvT[i,k] * diff(u[i,0], _p[k])
247 else:
248 raise RuntimeError("Taking divergence of a scalar.")
249 return d
250
251
252 -def curl(u, GinvT=None):
253 g = grad(u, GinvT)
254 c = swiginac.matrix(2,1)
255 c[0,0] = -g[1,0]
256 c[1,0] = g[0,0]
257 return c
258
259
261 return div(grad(u, GinvT), GinvT)
262
263
264
265 if __name__ == '__main__':
266 a = swiginac.matrix(3,1,[1,0,0])
267 b = swiginac.matrix(3,1,[0,1,0])
268 c = swiginac.matrix(3,1,[0,0,1])
269
270 print ""
271 print cross(a,b)
272 print cross(a,c)
273 print cross(c,b)
274 print cross(b,c)
275
276 import SyFi
277 SyFi.initSyFi(2)
278
279 x, y, z = symbols(["x", "y", "z"])
280 x2, y2, z2 = symbols(["x", "y", "z"])
281 assert (x-x2).is_zero()
282 assert (y-y2).is_zero()
283 assert (z-z2).is_zero()
284
285 p2 = swiginac.matrix(2, 1, [x,y])
286 p3 = swiginac.matrix(3, 1, [x,y,z])
287
288
289 A = swiginac.matrix(2,2, [x*x*x, x*x*y, x*y*y, y*y*y])
290 sh = shape(A)
291 assert rank(A) == 2
292 assert len(sh) == 2
293 assert sh[0] == 2
294 assert sh[1] == 2
295 print ""
296 print "A", A
297 print "div A", div(A)
298
299
300 v = swiginac.matrix(2, 1, [3*x,5*y])
301 sh = shape(v)
302 assert rank(v) == 1
303 assert len(sh) == 1
304 assert sh[0] == 2
305 print ""
306 print "v", v
307 print "div v", div(v)
308 print "grad v", grad(v)
309 print "dv/dx ", diff(v, p2)
310
311
312 v = swiginac.matrix(2, 1, [3*y,5*x])
313 sh = shape(v)
314 assert rank(v) == 1
315 assert len(sh) == 1
316 assert sh[0] == 2
317 print ""
318 print "v", v
319 print "div v", div(v)
320 print "grad v", grad(v)
321 print "dv/dx ", diff(v, p2)
322