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

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