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

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