Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /branches/ein16/synth/d2/test_eval.py
ViewVC logotype

Annotation of /branches/ein16/synth/d2/test_eval.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4237 - (view) (download) (as text)

1 : cchiw 3915 import sympy
2 :     from sympy import *
3 :     #symbols
4 :     x,y,z =symbols('x y z')
5 : cchiw 3865 import sys
6 :     import re
7 : cchiw 4188 import math
8 : cchiw 3915 from obj_ty import *
9 :     from obj_apply import *
10 :     from obj_operator import *
11 : cchiw 3939 from obj_field import *
12 : cchiw 3865
13 : cchiw 3874
14 : cchiw 3865
15 : cchiw 3915 # *************************** unary operators ***************************
16 : cchiw 3939 # binary operators
17 :     def fn_add(exp1,exp2):
18 :     return exp1+exp2
19 :     def fn_subtract(exp1,exp2):
20 :     return exp1-exp2
21 :     # scaling operator
22 : cchiw 3998 def fn_multiplication(exp_s, t):
23 : cchiw 3939 exp_t = field.get_data(t)
24 :     ityp1 = field.get_ty(t)
25 :     shape1 = fty.get_shape(ityp1)
26 : cchiw 3946 #print "exp_s",exp_s,"exp_t",exp_t
27 : cchiw 3939 if(field.is_Scalar(t)):
28 :     return exp_s* exp_t
29 :     elif(field.is_Vector(t)):
30 :     [n1] = shape1 #vector
31 :     rtn = []
32 :     for i in range(n1):
33 :     rtn.append(exp_s*exp_t[i])
34 :     return rtn
35 :     elif(field.is_Matrix(t)):
36 :     [n1,n2] = shape1
37 :     rtn = []
38 :     for i in range(n1):
39 : cchiw 3946 tmp = []
40 : cchiw 3939 for j in range(n2):
41 : cchiw 3946 tmp.append(exp_s*exp_t[i][j])
42 :     rtn.append(tmp)
43 : cchiw 3939 return rtn
44 :     elif(field.is_Ten3(t)):
45 :     [n1,n2,n3] = shape1
46 :     rtn = []
47 :     for i in range(n1):
48 : cchiw 3946 tmpI = []
49 : cchiw 3939 for j in range(n2):
50 : cchiw 3946 tmpJ = []
51 : cchiw 3939 for k in range(n3):
52 : cchiw 3946 tmpJ.append(exp_s*exp_t[i][j][k])
53 :     tmpI.append(tmpJ)
54 :     rtn.append(tmpI)
55 : cchiw 3939 return rtn
56 :     else:
57 :     raise "unsupported scaling"
58 : cchiw 3865
59 : cchiw 3998 #scaling of a field
60 :     def fn_scaling(fld1, fld2):
61 :     def get_sca():
62 :     if(field.is_Scalar(fld1)):
63 :     return (fld1, fld2)
64 :     else:
65 :     return (fld2, fld1)
66 :     (s, t) = get_sca()
67 :     exp_s = field.get_data(s)
68 :     return fn_multiplication(exp_s, t)
69 :    
70 :     #division of a field
71 :     def fn_division(t, s):
72 :     if(field.is_Scalar(s)):
73 :     exp_s = (1/field.get_data(s))
74 :     return fn_multiplication(exp_s, t)
75 :     else:
76 :     raise Exception ("err second arg in division should be a scalar")
77 :    
78 : cchiw 4236
79 :     # sine of field
80 :     def fn_negation(exp):
81 :     return -1*exp
82 :    
83 : cchiw 3946 # negation of field
84 :     def fn_negation(exp):
85 :     return -1*exp
86 :    
87 :    
88 : cchiw 3939 #evaluate cross product
89 :     def fn_cross(fld1, fld2):
90 :     exp1 = field.get_data(fld1)
91 :     ityp1 = field.get_ty(fld1)
92 :     exp2 = field.get_data(fld2)
93 :     ityp2 = field.get_ty(fld2)
94 :     #print " exp1: ",exp1," exp2: ",exp2
95 :     # vectors
96 :     n1 = fty.get_vecLength(ityp1) #vector
97 :     n2 = fty.get_vecLength(ityp2)
98 :     if(n1==2):
99 :     return (exp1[0]*exp2[1]) -(exp1[1]*exp2[0])
100 :     elif(n1==3):
101 :     x0= (exp1[1]*exp2[2]) -(exp1[2]*exp2[1])
102 :     x1= (exp1[2]*exp2[0]) -(exp1[0]*exp2[2])
103 :     x2= (exp1[0]*exp2[1]) -(exp1[1]*exp2[0])
104 :     return [x0,x1,x2]
105 :     else:
106 :     raise "unsupported type for cross product"
107 :    
108 : cchiw 3915 #gradient of field
109 :     def fn_grad(exp, dim):
110 :     if (dim==1):
111 : cchiw 4230 return diff(exp,x)
112 : cchiw 3915 elif (dim==2):
113 :     return [diff(exp,x), diff(exp,y)]
114 :     elif (dim==3):
115 :     return [diff(exp,x), diff(exp,y), diff(exp,z)]
116 :     else:
117 :     raise "dimension not supported"
118 : cchiw 3874
119 : cchiw 3939 #evaluate divergence
120 :     def fn_divergence(fld):
121 :     exp = field.get_data(fld)
122 :     ityp = field.get_ty(fld)
123 :     #print " exp1: ",exp1," exp2: ",exp2
124 :     # vectors
125 :     n1 = fty.get_vecLength(ityp) #vector
126 :     if(n1==2):
127 :     return diff(exp[0],x)+diff(exp[1],y)
128 :     elif(n1==3):
129 :    
130 :     return diff(exp[0],x)+diff(exp[1],y)+diff(exp[2],z)
131 :     else:
132 :     raise "unsupported type for divergence"
133 : cchiw 3874
134 : cchiw 3939 #evaluate cross product
135 :     def fn_curl(fld):
136 :     exp = field.get_data(fld)
137 :     ityp = field.get_ty(fld)
138 :     dim = field.get_dim(fld)
139 :     n = fty.get_vecLength(ityp) #vector
140 :     if(n!=dim):
141 :     raise (" type not supported for curl")
142 :     if(n==2):
143 :     return diff(exp[1], x) - diff(exp[0], y)
144 :     elif(n==3):
145 :     x0= diff(exp[2],y) - diff(exp[1],z)
146 :     x1= diff(exp[0],z) - diff(exp[2],x)
147 :     x2= diff(exp[1],x) - diff(exp[0],y)
148 :     return [x0,x1,x2]
149 :     else:
150 :     raise "unsupported type for cross product"
151 : cchiw 3874
152 : cchiw 3939 #evaluate jacob
153 :     def fn_jacob(fld):
154 :     exp = field.get_data(fld)
155 :     ityp = field.get_ty(fld)
156 :     dim = field.get_dim(fld)
157 :     #print " exp1: ",exp1," exp2: ",exp2
158 :     # vectors
159 :     n = fty.get_vecLength(ityp) #vector
160 :     if(n!=dim):
161 :     raise (" type not supported for jacob")
162 :     if(n==2):
163 : cchiw 3946 return [[diff(exp[0],x), diff(exp[0],y)],
164 :     [diff(exp[1],x), diff(exp[1],y)]]
165 : cchiw 3939 elif(n==3):
166 : cchiw 3946 return [[diff(exp[0],x), diff(exp[0],y), diff(exp[0],z)],
167 :     [diff(exp[1],x), diff(exp[1],y), diff(exp[1],z)],
168 :     [diff(exp[2],x), diff(exp[2],y), diff(exp[2],z)]]
169 : cchiw 3939 else:
170 : cchiw 3998 raise "unsupported type for jacob"
171 : cchiw 3939
172 : cchiw 3998 #evaluate norm
173 :     def fn_norm(fld, dim):
174 :     exp = field.get_data(fld)
175 :     ityp = field.get_ty(fld)
176 :     dim = field.get_dim(fld)
177 :     #print " exp1: ",exp1," exp2: ",exp2
178 :     # vectors
179 :     def iter (es):
180 :     sum = 0
181 :     for i in es:
182 :     sum+=i*i
183 : cchiw 4210 #print "\nsum",sum
184 : cchiw 4188 rtn = (sum)**0.5
185 : cchiw 4210 #print "\nrtn",rtn
186 : cchiw 4188 return rtn
187 : cchiw 3998 if(field.is_Scalar(fld)):
188 :     [] = fty.get_shape(ityp)
189 :     return exp
190 :     elif(field.is_Vector(fld)):
191 :     [n] = fty.get_shape(ityp)
192 :     rtn = []
193 :     for i in range(n):
194 :     rtn.append(exp[i])
195 :     return iter(rtn)
196 :     elif(field.is_Matrix(fld)):
197 :     [n, m] = fty.get_shape(ityp)
198 :     rtn = []
199 :     for i in range(n):
200 :     for j in range(m):
201 :     rtn.append(exp[i][j])
202 :     return iter(rtn)
203 :     elif(field.is_Ten3(fld)):
204 :     [n, m, p] = fty.get_shape(ityp)
205 :     rtn = []
206 :     for i in range(n):
207 :     for j in range(m):
208 :     for k in range(p):
209 :     rtn.append(exp[i][j][k])
210 :     return iter(rtn)
211 :     else:
212 :     raise "unsupported type for norm"
213 : cchiw 3939
214 : cchiw 3998 #evaluate norm
215 :     def fn_normalize(fld, dim):
216 :     exp = field.get_data(fld)
217 :     ityp = field.get_ty(fld)
218 :     dim = field.get_dim(fld)
219 :     #print " exp1: ",exp1," exp2: ",exp2
220 :     norm = fn_norm(fld, dim)
221 :     if(field.is_Scalar(fld)):
222 :     #print "scal",exp
223 :     return exp
224 :     elif(field.is_Vector(fld)):
225 :     [n] = fty.get_shape(ityp)
226 :     rtn = []
227 :     for i in range(n):
228 :     rtn.append(exp[i]/norm)
229 :     #print "vec",rtn
230 :     return rtn
231 :     elif(field.is_Matrix(fld)):
232 :     [n, m] = fty.get_shape(ityp)
233 :     rtn = []
234 :     for i in range(n):
235 :     rtni = []
236 :     for j in range(m):
237 :     rtni.append(exp[i][j]/norm)
238 :     rtn.append(rtni)
239 :     #print "matrix:",rtn
240 :     return rtn
241 :     elif(field.is_Ten3(fld)):
242 :     [n, m, p] = fty.get_shape(ityp)
243 :     rtn = []
244 :     for i in range(n):
245 :     rtni = []
246 :     for j in range(m):
247 :     rtnj = []
248 :     for k in range(p):
249 :     rtnj.append(exp[i][j][k]/norm)
250 :     rtni.append( rtnj)
251 :     rtn.append( rtni)
252 :     #print "ten3",rtn
253 :     return rtn
254 :     else:
255 :     raise "unsupported type for norm"
256 :    
257 : cchiw 4237 #[0]
258 :     def fn_slicev0(fld1):
259 :     exp1 = field.get_data(fld1)
260 :     ityp1 = field.get_ty(fld1)
261 :     return exp1[0]
262 :    
263 :    
264 : cchiw 3998 #evaluate slice
265 : cchiw 4237 #[1]
266 :     def fn_slicev1(fld1):
267 : cchiw 3998 exp1 = field.get_data(fld1)
268 :     ityp1 = field.get_ty(fld1)
269 : cchiw 4237 return exp1[1]
270 :    
271 :     #evaluate slice
272 :     #[1,:]
273 :     def fn_slicem0(fld1):
274 :     exp1 = field.get_data(fld1)
275 :     ityp1 = field.get_ty(fld1)
276 : cchiw 3998 rtn=[]
277 :     if(fty.is_Matrix(ityp1)):
278 :     [n2,n3] = fty.get_shape(ityp1)
279 : cchiw 4237 for i in range(n3):
280 :     rtn.append(exp1[1][i])
281 : cchiw 3998 return rtn
282 :     else:
283 :     raise "unsupported type for slice"
284 :    
285 : cchiw 4237 #evaluate slice
286 :     #[:,0]
287 :     def fn_slicem1(fld1):
288 :     exp1 = field.get_data(fld1)
289 :     ityp1 = field.get_ty(fld1)
290 :     rtn=[]
291 :     if(fty.is_Matrix(ityp1)):
292 :     [n2,n3] = fty.get_shape(ityp1)
293 :     for i in range(n2):
294 :     rtn.append(exp1[i][0])
295 :     return rtn
296 :     else:
297 :     raise "unsupported type for slice"
298 :    
299 :    
300 : cchiw 3998 #evaluate trace
301 :     def fn_trace(fld):
302 :     exp = field.get_data(fld)
303 :     ityp = field.get_ty(fld)
304 :     rtn=[]
305 :     if(field.is_Matrix(fld)):
306 :     [n, m] = fty.get_shape(ityp)
307 :     if (n!=m):
308 :     raise Exception("matrix is not identitical")
309 :     rtn = exp[0][0]+exp[1][1]
310 :     if(n==2):
311 :     return rtn
312 :     elif(n==3):
313 :     return rtn+exp[2][2]
314 :     else:
315 :     raise "unsupported matrix field"
316 :     else:
317 :     raise "unsupported trace"
318 :    
319 :     #evaluate transpose
320 :     def fn_transpose(fld):
321 :     exp = field.get_data(fld)
322 :     ityp = field.get_ty(fld)
323 :     if(field.is_Matrix(fld)):
324 :     [n, m] = fty.get_shape(ityp)
325 :     rtn = []
326 :     for i in range(n):
327 :     rtni = []
328 :     for j in range(m):
329 :     rtni.append(exp[j][i])
330 :     rtn.append(rtni)
331 :     return rtn
332 :     else:
333 :     raise "unsupported transpose"
334 :    
335 :     #evaluate det
336 :     def fn_det(fld):
337 :     exp = field.get_data(fld)
338 :     ityp = field.get_ty(fld)
339 :     rtn=[]
340 :     if(field.is_Matrix(fld)):
341 :     [n, m] = fty.get_shape(ityp)
342 :     if (n!=m):
343 :     raise Exception("matrix is not identitical")
344 :     a = exp[0][0]
345 :     d = exp[1][1]
346 :     c = exp[1][0]
347 :     b = exp[0][1]
348 :     if(n==2):
349 :     return a*d-b*c
350 :     elif(n==3):
351 :     a = exp[0][0]
352 :     b = exp[0][1]
353 :     c = exp[0][2]
354 :     d = exp[1][0]
355 :     e = exp[1][1]
356 :     f = exp[1][2]
357 :     g = exp[2][0]
358 :     h = exp[2][1]
359 :     i = exp[2][2]
360 :     return a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)
361 :     else:
362 :     raise "unsupported matrix field"
363 :     else:
364 :     raise "unsupported trace"
365 :    
366 :    
367 : cchiw 3939 #evaluate outer product
368 :     def fn_outer(fld1, fld2):
369 :     exp1 = field.get_data(fld1)
370 :     ityp1 = field.get_ty(fld1)
371 :     exp2 = field.get_data(fld2)
372 :     ityp2 = field.get_ty(fld2)
373 :     rtn=[]
374 : cchiw 3946 #print "exp1",exp1,"ityp1",ityp1.name,"-length",len(exp1)
375 :     #print "exp2",exp2,"ityp2",ityp2.name,"-length",len(exp2)
376 :    
377 : cchiw 3939 if(fty.is_Vector(ityp1)):
378 :     n1= fty.get_vecLength(ityp1)
379 :     if(fty.is_Vector(ityp2)):
380 :     #both vectors
381 :     n2= fty.get_vecLength(ityp2)
382 :     for i in range(n1):
383 : cchiw 3946 tmpI = []
384 : cchiw 3939 for j in range(n2):
385 : cchiw 3946 tmpI.append(exp1[i]*exp2[j])
386 :     rtn.append(tmpI)
387 : cchiw 3939 return rtn
388 :     elif(fty.is_Matrix(ityp2)):
389 :     [n2,n3] = fty.get_shape(ityp2)
390 :     for i in range(n1):
391 : cchiw 3946 tmpI = []
392 : cchiw 3939 for j in range(n2):
393 : cchiw 3946 tmpJ = []
394 : cchiw 3939 for k in range(n3):
395 : cchiw 3946 tmpJ.append(exp1[i]*exp2[j][k])
396 :     tmpI.append(tmpJ)
397 :     rtn.append(tmpI)
398 : cchiw 3939 return rtn
399 :     elif(fty.is_Matrix(ityp1)):
400 :     [n1,n2] = fty.get_shape(ityp1)
401 :     if(fty.is_Vector(ityp2)):
402 :     n3= fty.get_vecLength(ityp2)
403 :     for i in range(n1):
404 : cchiw 3946 tmpI = []
405 : cchiw 3939 for j in range(n2):
406 : cchiw 3946 tmpJ = []
407 : cchiw 3939 for k in range(n3):
408 : cchiw 3946 tmpJ.append(exp1[i][j]*exp2[k])
409 :     tmpI.append(tmpJ)
410 :     rtn.append(tmpI)
411 : cchiw 3939 return rtn
412 :     elif(fty.is_Matrix(ityp2)):
413 :     [n3,n4] = fty.get_shape(ityp2)
414 :     for i in range(n1):
415 : cchiw 3946 tmpI = []
416 : cchiw 3939 for j in range(n2):
417 : cchiw 3946 tmpJ = []
418 : cchiw 3939 for k in range(n3):
419 : cchiw 3946 tmpK = []
420 : cchiw 3939 for l in range(n4):
421 : cchiw 3946 tmpK.append(exp1[i][j]*exp2[k][l])
422 :     tmpJ.append(tmpK)
423 :     tmpI.append(tmpJ)
424 :     rtn.append(tmpI)
425 : cchiw 3939 return rtn
426 :     else:
427 :     raise "outer product is not supported"
428 :     else:
429 :     raise "outer product is not supported"
430 :     #evaluate inner product
431 :     def fn_inner(fld1, fld2):
432 :     exp1 = field.get_data(fld1)
433 :     ityp1 = field.get_ty(fld1)
434 :     exp2 = field.get_data(fld2)
435 :     ityp2 = field.get_ty(fld2)
436 :     #print " exp1: ",exp1," exp2: ",exp2
437 :     # vectors
438 :     if(fty.is_Vector(ityp1)):
439 :     n1 = fty.get_vecLength(ityp1) #vector
440 :     if(fty.is_Vector(ityp2)):
441 :     #length of vetors
442 :     rtn=0
443 :     n2 = fty.get_vecLength(ityp2)
444 :     for s in range(n1):
445 :     curr = exp1[s]*exp2[s]
446 :     #print (" exp1[s]: ",exp1[s]," exp2[s]: ",exp2[s],"cur",curr)
447 :     rtn += curr
448 :     return rtn
449 :     elif(fty.is_Matrix(ityp2)):
450 :     [n2] = fty.drop_last(ityp2) #matrix
451 :     rtn=[]
452 :     for i in range(n2):
453 :     sumrtn=0
454 :     for s in range(n1):
455 :     curr = exp1[s]*exp2[s][i]
456 :     sumrtn += curr
457 :     rtn.append(sumrtn)
458 :     return rtn
459 :     elif(fty.is_Ten3(ityp2)):
460 :     [n2,n3] = fty.drop_last(ityp2)
461 : cchiw 3946 rtn = []
462 : cchiw 3939 for i in range(n2):
463 : cchiw 3946 tmpJ = []
464 : cchiw 3939 for j in range(n3):
465 :     sumrtn=0
466 :     for s in range(n1):
467 :     curr = exp1[s]*exp2[s][i][j]
468 :     sumrtn += curr
469 : cchiw 3946 tmpJ.append(sumrtn)
470 :     rtn.append(tmpJ)
471 :     return rtn
472 : cchiw 3939 else:
473 :     raise "inner product is not supported"
474 :     elif(fty.is_Matrix(ityp1)):
475 :     n2 = fty.get_first_ix(ityp1) #matrix
476 :     if(fty.is_Vector(ityp2)):
477 :     ns = fty.get_vecLength(ityp2) #vector
478 :     rtn=[]
479 :     for i in range(n2):
480 :     sumrtn=0
481 :     for s in range(ns):
482 :     curr = exp1[i][s]*exp2[s]
483 :     sumrtn += curr
484 :     rtn.append(sumrtn)
485 :     return rtn
486 :     else:
487 :     raise "inner product is not supported"
488 :     elif(fty.is_Ten3(ityp1)):
489 :     [n1,n2] = fty.drop_first(ityp1)
490 :     if(fty.is_Vector(ityp2)):
491 :     ns = fty.get_vecLength(ityp2)
492 :     rtn=[]
493 :     for i in range(n1):
494 : cchiw 3946 tmpI=[]
495 : cchiw 3939 for j in range(n2):
496 :     sumrtn=0
497 :     for s in range(ns):
498 :     curr = exp1[i][j][s]*exp2[s]
499 :     sumrtn += curr
500 : cchiw 3946 tmpI.append(sumrtn)
501 :     rtn.append(tmpI)
502 : cchiw 3939 return rtn
503 :     else:
504 :     raise "inner product is not supported"
505 :     else:
506 :     raise "inner product is not supported"
507 :    
508 :    
509 : cchiw 3915 # *************************** generic apply operators ***************************
510 :     #unary operator on a vector
511 : cchiw 3939 def applyToVector(vec, unary):
512 : cchiw 3915 rtn = []
513 :     for v in vec:
514 :     rtn.append(unary(v))
515 :     return rtn
516 :     #binary operator on a vector
517 : cchiw 3939 def applyToVectors(vecA, vecB, binary):
518 : cchiw 3915 rtn = []
519 :     for (a,b) in zip(vecA,vecB):
520 :     x= binary(a,b)
521 :     rtn.append(x)
522 :     return rtn
523 : cchiw 3874
524 : cchiw 3946 def applyToM(vec, unary):
525 :     rtn = []
526 :     for i in vec:
527 :     tmpI = []
528 :     for v in i:
529 :     tmpI.append(unary(v))
530 :     rtn.append(tmpI)
531 :     return rtn
532 :    
533 : cchiw 4230 def applyToMs(vecA,vecB, unary):
534 :     rtn = []
535 :     for (a,b) in zip(vecA,vecB):
536 :     tmpI = []
537 :     for (u,v) in zip(a,b):
538 :     tmpI.append(unary(u, v))
539 :     rtn.append(tmpI)
540 :     return rtn
541 :    
542 :    
543 : cchiw 3946 def applyToT3(vec, unary):
544 :     rtn = []
545 :     for i in vec:
546 :     tmpI = []
547 :     for j in i:
548 :     tmpJ = []
549 :     for v in j:
550 :     tmpJ.append(unary(v))
551 :     tmpI.append(tmpJ)
552 :     rtn.append(tmpI)
553 :     return rtn
554 :    
555 : cchiw 3915 # *************************** apply to scalars or vectors ***************************
556 :     #apply operators to expression
557 :     # return output types and expression
558 :     # unary operator
559 :     # exp: scalar types
560 : cchiw 3874
561 : cchiw 3939 def applyToExp_U_S(fn_name, fld):
562 :     exp = field.get_data(fld)
563 :     dim = field.get_dim(fld)
564 :     #print fn_name
565 : cchiw 3915 if(op_probe==fn_name): #probing
566 : cchiw 3939 return exp
567 : cchiw 3915 elif(op_negation==fn_name): #negation
568 : cchiw 3939 return fn_negation(exp)
569 : cchiw 3915 elif(op_gradient==fn_name): #gradient
570 : cchiw 3939 return fn_grad(exp, dim)
571 : cchiw 3915 else:
572 : cchiw 3939 raise Exception("unsupported unary operator on scalar field:"+ fn_name.name)
573 : cchiw 3874
574 : cchiw 3939 # unary operator
575 :     # exp: vector types
576 :     def applyToExp_U_V(fn_name, fld):
577 :     exp = field.get_data(fld)
578 :     if(op_probe==fn_name): #probing
579 :     return exp
580 :     elif(op_negation==fn_name): #negation
581 :     return applyToVector(exp, fn_negation)
582 :     elif(op_divergence==fn_name):
583 :     return fn_divergence(fld)
584 :     elif(op_curl==fn_name):
585 :     return fn_curl(fld)
586 :     elif(op_jacob==fn_name): #jacob
587 :     return fn_jacob(fld)
588 : cchiw 4237 elif(op_slicev0==fn_name) :
589 :     return fn_slicev0(fld)
590 :     elif(op_slicev1==fn_name):
591 :     return fn_slicev1(fld)
592 : cchiw 3939 else:
593 :     raise Exception("unsupported unary operator:"+ fn_name.name)
594 : cchiw 3874
595 : cchiw 3946 def applyToExp_U_M(fn_name, fld):
596 :     exp = field.get_data(fld)
597 :     if(op_probe==fn_name): #probing
598 :     return exp
599 :     elif(op_negation==fn_name): #negation
600 :     return applyToM(exp, fn_negation)
601 :     elif(op_jacob==fn_name): #jacob
602 :     return fn_jacob(fld)
603 : cchiw 4237 elif(op_slicem0==fn_name) :
604 :     return fn_slicem0(fld)
605 :     elif(op_slicem1==fn_name):
606 :     return fn_slicem1(fld)
607 : cchiw 3998 elif(op_trace == fn_name):
608 :     return fn_trace(fld)
609 :     elif(op_transpose==fn_name):
610 :     return fn_transpose(fld)
611 :     elif(op_det==fn_name):
612 :     return fn_det(fld)
613 : cchiw 3946 else:
614 :     raise Exception("unsupported unary operator:"+ fn_name.name)
615 :    
616 : cchiw 3939 def applyToExp_U_T3(fn_name, fld):
617 :     exp = field.get_data(fld)
618 :     if(op_probe==fn_name): #probing
619 :     return exp
620 :     elif(op_negation==fn_name): #negation
621 : cchiw 3946
622 :     return applyToT3(exp, fn_negation)
623 : cchiw 3939 elif(op_jacob==fn_name): #jacob
624 :     return fn_jacob(fld)
625 :     else:
626 :     raise Exception("unsupported unary operator:"+ fn_name.name)
627 :    
628 : cchiw 3915 # binary operator
629 :     # exp: scalar types
630 :     def applyToExp_B_S(e):
631 :     fn_name=e.opr
632 : cchiw 3939 (fld1,fld2) = apply.get_binary(e)
633 :     exp1 = field.get_data(fld1)
634 :     exp2 = field.get_data(fld2)
635 :     #print fn_name
636 : cchiw 3915 if(op_add==fn_name):#addition
637 : cchiw 3939 return fn_add(exp1,exp2)
638 : cchiw 3915 elif(op_subtract==fn_name):#subtract
639 : cchiw 3939 return fn_subtract(exp1,exp2)
640 :     elif(op_scale==fn_name): #scaling
641 :     return fn_scaling(fld1,fld2)
642 : cchiw 3998 elif(op_division==fn_name): #division
643 :     return fn_division(fld1,fld2)
644 : cchiw 3915 else:
645 : cchiw 3939 raise Exception("unsupported binary operator on scalar fields:"+ fn_name.name)
646 : cchiw 3874
647 : cchiw 3939 # binary, args do not need to have the same shape
648 :     def applyToExp_B_uneven(e):
649 : cchiw 3915 fn_name=e.opr
650 : cchiw 3939 (fld1,fld2) = apply.get_binary(e)
651 :     exp1 = field.get_data(fld1)
652 :     exp2 = field.get_data(fld2)
653 :     if(op_outer==fn_name):
654 :     return fn_outer(fld1, fld2)
655 :     elif(op_inner==fn_name):
656 :     return fn_inner(fld1, fld2)
657 :     elif(op_scale==fn_name): #scaling
658 :     return fn_scaling(fld1,fld2)
659 : cchiw 3915 else:
660 :     raise Exception("unsupported unary operator:",op_name)
661 : cchiw 3865
662 : cchiw 3915
663 :     # binary operator
664 : cchiw 3939 # args have the same shape
665 : cchiw 3915 def applyToExp_B_V(e):
666 :     fn_name=e.opr
667 : cchiw 3939 (fld1,fld2) = apply.get_binary(e)
668 :     exp1 = field.get_data(fld1)
669 :     exp2 = field.get_data(fld2)
670 : cchiw 3915 if(op_add==fn_name):#addition
671 : cchiw 3939 return applyToVectors(exp1, exp2, fn_add)
672 : cchiw 3915 elif(op_subtract==fn_name):#subtract
673 : cchiw 3939 return applyToVectors(exp1, exp2, fn_subtract)
674 :     elif(op_cross==fn_name):
675 :     return fn_cross(fld1, fld2)
676 : cchiw 3915 else:
677 : cchiw 3939 return applyToExp_B_uneven(e)
678 : cchiw 3915
679 : cchiw 4230 def applyToExp_B_M(e):
680 :     fn_name=e.opr
681 :     (fld1,fld2) = apply.get_binary(e)
682 :     exp1 = field.get_data(fld1)
683 :     exp2 = field.get_data(fld2)
684 :     if(op_add==fn_name):#addition
685 :     return applyToMs(exp1, exp2, fn_add)
686 :     elif(op_subtract==fn_name):#subtract
687 :     return applyToMs(exp1, exp2, fn_subtract)
688 :     else:
689 :     return applyToExp_B_uneven(e)
690 : cchiw 3939
691 : cchiw 4230
692 : cchiw 3915 # *************************** unary/binary operators ***************************
693 :     def unary(e):
694 : cchiw 3939 #apply.toStr(e)
695 :     fld =apply.get_unary(e)
696 :     fn_name=e.opr
697 :     exp = field.get_data(fld)
698 :     dim = field.get_dim(fld)
699 : cchiw 3998 if(op_norm==fn_name):#norm
700 :     return fn_norm(fld, dim)
701 :     if(op_normalize==fn_name):#normalize
702 :     x= fn_normalize(fld, dim)
703 :     return x
704 :     elif (field.is_Scalar(fld)): # input is a scalar field
705 : cchiw 3939 return applyToExp_U_S(fn_name, fld)
706 :     elif(field.is_Vector(fld)): # input is a vector field
707 :     return applyToExp_U_V(fn_name, fld)
708 : cchiw 3946 elif(field.is_Matrix(fld)): # input is a vector field
709 :     return applyToExp_U_M(fn_name, fld)
710 : cchiw 3915 else:
711 : cchiw 3939 return applyToExp_U_T3(fn_name, fld)
712 : cchiw 3915
713 :     def binary(e):
714 : cchiw 3939 (f, g) =apply.get_binary(e)
715 : cchiw 3998 fn_name = e.opr
716 :     # type is checked elsewhere or does not matter
717 :     if(op_division==fn_name): #division
718 :     return fn_division(f, g)
719 :     elif (field.is_Scalar(f) and field.is_Scalar(g)): # input is a scalar field
720 : cchiw 3915 return applyToExp_B_S(e)
721 : cchiw 3939 elif (field.is_Vector(f)):# input is a vector field
722 :     if(field.is_Vector(g)):
723 :     return applyToExp_B_V(e)
724 :     else: # input is a vector field, _
725 : cchiw 4230 return applyToExp_B_V(e)
726 :     elif (field.is_Matrix(f)):# input is a matrix field
727 :     if(field.is_Matrix(g)):
728 :     return applyToExp_B_M(e)
729 :     else:
730 :     return applyToExp_B_V(e)
731 : cchiw 3915 else:
732 : cchiw 4230 return applyToExp_B_V(e)
733 : cchiw 3915
734 : cchiw 3939 def applyUnaryOnce(oexp_inner,app_inner,app_outer):
735 : cchiw 3998 #print "applyUnaryOnce"
736 : cchiw 3939 #apply.toStr(app_inner)
737 :     oty_inner = apply.get_oty(app_inner)
738 :     oty_outer = apply.get_oty(app_outer)
739 :     opr_outer = app_outer.opr
740 : cchiw 3946 #print "oexp_inner",oexp_inner,"opr_outer",opr_outer.name
741 : cchiw 3939 lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
742 : cchiw 3946 app_tmp = apply("tmp", opr_outer, lhs_tmp, None, oty_outer, true, true)
743 : cchiw 3939 oexp_tmp =unary(app_tmp)
744 : cchiw 3946 #print " oexp_tmp", oexp_tmp
745 : cchiw 3939 return (oty_outer, oexp_tmp)
746 :    
747 : cchiw 3946 def applyBinaryOnce(oexp_inner,app_inner,app_outer,rhs):
748 :     oty_inner = apply.get_oty(app_inner)
749 :     oty_outer = apply.get_oty(app_outer)
750 :     opr_outer = app_outer.opr
751 :    
752 :     lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
753 :    
754 :     app_tmp = apply("tmp", opr_outer, lhs_tmp, rhs, oty_outer, true,true)
755 :     oexp_tmp =binary(app_tmp)
756 :     return (oty_outer, oexp_tmp)
757 : cchiw 3939
758 : cchiw 3946
759 : cchiw 3915 # operators with scalar field and vector field
760 :     def sort(e):
761 : cchiw 3939 #apply.toStr(e)
762 :     arity = apply.get_arity(e)
763 : cchiw 3946 if(e.isrootlhs): # is root
764 : cchiw 3939 #print "sort is a root"
765 :     oty = apply.get_oty(e)
766 :     if (arity ==1):
767 :     return (oty, unary(e))
768 :     elif (arity ==2): # binary operator
769 :     return (oty, binary(e))
770 :     else:
771 :     raise Exception ("arity is not supported: "+str(arity_outer))
772 : cchiw 3915 else:
773 : cchiw 3939 app_outer = e
774 :     arity_outer = arity
775 : cchiw 3946 #print "app_outer",app_outer.opr.name
776 : cchiw 3939 if (arity_outer==1): #assuming both arity
777 :     app_inner = apply.get_unary(app_outer)
778 : cchiw 3946 #print "outer(1) app_inner:",app_inner.opr.name
779 : cchiw 3939 arity_inner= app_inner.opr.arity
780 :     if (arity_inner==1):
781 :     oexp_inner = unary(app_inner)
782 : cchiw 3946
783 : cchiw 3939 (oty_outer, oexp_tmp) = applyUnaryOnce(oexp_inner ,app_inner, app_outer)
784 : cchiw 3946
785 : cchiw 3939 return (oty_outer, oexp_tmp)
786 :     elif(arity_inner==2):
787 :     oexp_inner = binary(app_inner)
788 :     (oty_outer, oexp_tmp) = applyUnaryOnce(oexp_inner ,app_inner, app_outer)
789 : cchiw 3946
790 : cchiw 3939 return (oty_outer, oexp_tmp)
791 :     else:
792 :     raise Exception ("arity is not supported: "+str(arity_outer))
793 : cchiw 3946 elif (arity_outer==2):
794 :     (app_inner, G) = apply.get_binary(app_outer)
795 :     arity_inner= app_inner.opr.arity
796 :     #print "outer(2) app_inner",app_inner.opr.name
797 :     if (arity_inner==1):
798 :     oexp_inner = unary(app_inner)
799 :     rhs = G
800 :     (oty_outer, oexp_tmp) = applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
801 :     return (oty_outer, oexp_tmp)
802 :     elif(arity_inner==2):
803 :     oexp_inner = binary(app_inner)
804 :     #print "applying binary frst time", oexp_inner
805 :     rhs = G
806 :     (oty_outer, oexp_tmp) = applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
807 :     #print "applying binary second time", oexp_tmp
808 :     return (oty_outer, oexp_tmp)
809 :     else:
810 :     raise Exception ("arity is not supported: "+str(arity_outer))
811 : cchiw 3939 else:
812 :     raise Exception ("arity is not supported: "+str(arity_outer))
813 : cchiw 3915
814 :     # *************************** evaluate at positions ***************************
815 :     #evaluate scalar field exp
816 : cchiw 4158 def eval_d1(pos0, exp):
817 : cchiw 4230 #print "eval vec d1"
818 : cchiw 4158 #print "exp:",exp
819 : cchiw 4230 #print "pos0",pos0
820 : cchiw 4158 #evaluate field defined by coefficients at position
821 :     exp = exp.subs(x,pos0)
822 : cchiw 4230 #print "exp",exp
823 : cchiw 4158 return exp
824 :    
825 : cchiw 3915 def eval_d2(pos0, pos1, exp):
826 : cchiw 3998 #print "exp:",exp
827 : cchiw 3915 #evaluate field defined by coefficients at position
828 : cchiw 4158 exp = exp.subs(x,pos0)
829 :     exp = exp.subs(y,pos1)
830 :     return exp
831 : cchiw 3915
832 : cchiw 3939 def eval_d3(pos0, pos1, pos2, exp):
833 :     #evaluate field defined by coefficients at position
834 : cchiw 4158 exp = exp.subs(x,pos0)
835 :     exp = exp.subs(y,pos1)
836 :     exp = exp.subs(z,pos2)
837 :     return exp
838 : cchiw 3939
839 : cchiw 4158 #evaluate vector field [exp]
840 :     def eval_vec_d1(pos0, vec):
841 : cchiw 4230 #print "eval vec d1"
842 : cchiw 4158 rtn = []
843 :     for v in vec:
844 :     rtn.append(eval_d1(pos0, v))
845 :     return rtn
846 :    
847 : cchiw 3915 #evaluate vector field [exp,exp]
848 :     def eval_vec_d2(pos0, pos1, vec):
849 : cchiw 3939 #print "eval_vec_d2 vec:",vec
850 : cchiw 3915 rtn = []
851 :     for v in vec:
852 :     rtn.append(eval_d2(pos0, pos1, v))
853 :     return rtn
854 :    
855 : cchiw 4158 def eval_ten3_d1(pos0,ten3):
856 :     rtn = []
857 :     for i in ten3:
858 :     for j in i:
859 :     for v in j:
860 :     rtn.append(eval_d1(pos0, v))
861 :     return rtn
862 :    
863 :    
864 :    
865 : cchiw 3939 #evaluate vector field [exp,exp]
866 : cchiw 4158 def eval_mat_d1(pos0, mat):
867 :     #print "eval_vec_d2 vec:",vec
868 :     rtn = []
869 :     for i in mat:
870 :     for v in i:
871 :     rtn.append(eval_d1(pos0, v))
872 :     return rtn
873 :    
874 :     #evaluate vector field [exp,exp]
875 : cchiw 3946 def eval_mat_d2(pos0, pos1, mat):
876 :     #print "eval_vec_d2 vec:",vec
877 :     rtn = []
878 : cchiw 4210 #print "eval_mat_d2 mat",mat
879 : cchiw 3946 for i in mat:
880 :     for v in i:
881 :     rtn.append(eval_d2(pos0, pos1, v))
882 :     return rtn
883 :    
884 :     def eval_ten3_d2(pos0, pos1, ten3):
885 :     #print "eval_vec_d2 vec:",vec
886 :     rtn = []
887 :     for i in ten3:
888 :     for j in i:
889 :     for v in j:
890 :     rtn.append(eval_d2(pos0, pos1, v))
891 :     return rtn
892 :    
893 :    
894 :    
895 :     #evaluate vector field [exp,exp]
896 : cchiw 3939 def eval_vec_d3(pos0, pos1, pos2, vec):
897 :     rtn = []
898 :     for v in vec:
899 :     rtn.append(eval_d3(pos0, pos1, pos2, v))
900 :     return rtn
901 :    
902 :    
903 : cchiw 3946 #evaluate vector field [exp,exp]
904 :     def eval_mat_d3(pos0, pos1, pos2, mat):
905 :     #print "eval_vec_d2 vec:",vec
906 :     rtn = []
907 :     for i in mat:
908 :     for v in i:
909 :     rtn.append(eval_d3(pos0, pos1, pos2, v))
910 :     return rtn
911 :    
912 :     def eval_ten3_d3(pos0, pos1, pos2,ten3):
913 :     rtn = []
914 :     for i in ten3:
915 :     for j in i:
916 :     for v in j:
917 :     rtn.append(eval_d3(pos0, pos1, pos2, v))
918 :     return rtn
919 :    
920 :    
921 : cchiw 4158
922 :     def iter_d1(k, pos, exp):
923 :     corr = []
924 :     for x in pos:
925 :     val = k(x, exp)
926 :     corr.append(val)
927 :     return corr
928 :    
929 : cchiw 3939 def iter_d2(k, pos, exp):
930 : cchiw 3874 corr = []
931 : cchiw 3939 #print "iter expr:", exp
932 :     #print "pos", pos
933 : cchiw 3874 for p in pos:
934 : cchiw 3939 #print "p", p
935 : cchiw 3874 x=p[0]
936 :     y=p[1]
937 : cchiw 3915 val = k(x,y,exp)
938 : cchiw 3874 corr.append(val)
939 :     return corr
940 : cchiw 3865
941 : cchiw 3939 def iter_d3(k, pos, exp):
942 :     corr = []
943 :     #print "iter exp:", exp
944 :     for p in pos:
945 :     x=p[0]
946 :     y=p[1]
947 :     z=p[2]
948 :     val = k(x,y,z, exp)
949 :     #print "pos: ",x,y,z, " val:", val
950 :     corr.append(val)
951 :     return corr
952 :    
953 : cchiw 3915 def probeField(otyp1,pos, ortn):
954 : cchiw 3939 dim = fty.get_dim(otyp1)
955 : cchiw 3946 #print "output type"+otyp1.name
956 : cchiw 4158 if (dim==1):
957 : cchiw 3939 def get_k():
958 :     if (fty.is_ScalarField(otyp1)): # output is a scalar field
959 : cchiw 4230 #print "s_d1"
960 : cchiw 4158 return eval_d1
961 :     elif (fty.is_VectorField(otyp1)):
962 : cchiw 4230 #print "v_d1"
963 : cchiw 4158 return eval_vec_d1
964 :     elif (fty.is_Matrix(otyp1)):
965 :     return eval_mat_d1
966 :     elif(fty.is_Ten3(otyp1)):
967 :     return eval_ten3_d1
968 :     else:
969 :     raise "error"+otyp1.name
970 :     return iter_d1(get_k(), pos, ortn)
971 :     elif (dim==2):
972 :     def get_k():
973 :     if (fty.is_ScalarField(otyp1)): # output is a scalar field
974 : cchiw 3939 return eval_d2
975 : cchiw 3946 elif (fty.is_VectorField(otyp1)):
976 :     return eval_vec_d2
977 :     elif (fty.is_Matrix(otyp1)):
978 :     return eval_mat_d2
979 :     elif(fty.is_Ten3(otyp1)):
980 :     return eval_ten3_d2
981 : cchiw 3939 else:
982 : cchiw 3946 raise "error"+otyp1.name
983 : cchiw 3939 return iter_d2(get_k(), pos, ortn)
984 :     elif (dim==3):
985 :     def get_k():
986 :     if (fty.is_ScalarField(otyp1)): # output is a scalar field
987 :     return eval_d3
988 : cchiw 3946 elif (fty.is_VectorField(otyp1)):
989 :     return eval_vec_d3
990 :     elif (fty.is_Matrix(otyp1)):
991 :     return eval_mat_d3
992 :     elif(fty.is_Ten3(otyp1)):
993 :     return eval_ten3_d3
994 : cchiw 3939 else:
995 : cchiw 3946 raise "error"+otyp1.name
996 : cchiw 3939 return iter_d3(get_k(), pos, ortn)
997 : cchiw 3865 else:
998 : cchiw 3939 raise "unsupported dimension"
999 : cchiw 3865
1000 : cchiw 3915 # *************************** main ***************************
1001 :    
1002 :     # operators with scalar field and vector field
1003 :     def eval(app, pos):
1004 : cchiw 3946 #print "evalname",app.name
1005 : cchiw 3939 #apply.toStr(app)
1006 : cchiw 3915 (otyp1, ortn) = sort(app) #apply operations to expressions
1007 : cchiw 4210 #print "ortn",ortn
1008 : cchiw 4188 rtn = probeField(otyp1, pos, ortn) #evaluate expression at positions
1009 : cchiw 4210 #print "rtn", rtn
1010 : cchiw 4188 return rtn

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0