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 4456 - (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 4456 from write import *
13 : cchiw 3865
14 : cchiw 3874
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 3939 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 : cchiw 4230 return diff(exp,x)
108 : cchiw 3915 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 : cchiw 4411 #print "inside jacob"
151 : cchiw 3939 exp = field.get_data(fld)
152 :     ityp = field.get_ty(fld)
153 :     dim = field.get_dim(fld)
154 :     #print " exp1: ",exp1," exp2: ",exp2
155 : cchiw 4411 shape = fty.get_shape(ityp)
156 : cchiw 3939 # vectors
157 : cchiw 4411 if(field.is_VectorField (fld)):
158 :     print "found a vector"
159 :     [n] = shape #vector
160 :     if(n!=dim):
161 :     raise (" type not supported for jacob")
162 :     if(dim==2):
163 :     return [[diff(exp[0],x), diff(exp[0],y)],
164 :     [diff(exp[1],x), diff(exp[1],y)]]
165 :     elif(dim==3):
166 :     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 :     else:
170 :     raise (" type not supported for jacob")
171 :     elif(field.is_Matrix (fld)):
172 :     #print "inside matrix"
173 :     [n,m] = shape #vector
174 :     if(n!=dim and n!=m):
175 :     raise (" type not supported for jacob")
176 :     elif(dim==2):
177 :     e_00 = exp[0][0]
178 :     e_10 = exp[1][0]
179 :     e_01 = exp[0][1]
180 :     e_11 = exp[1][1]
181 :     f_00 = [diff(e_00,x), diff(e_00,y)]
182 :     f_01 = [diff(e_01,x), diff(e_01,y)]
183 :     f_10 = [diff(e_10,x), diff(e_10,y)]
184 :     f_11 = [diff(e_11,x), diff(e_11,y)]
185 :     return [[f_00,f_01],[f_10,f_11]]
186 :     elif(dim==3):
187 :     e_00 = exp[0][0]
188 :     e_10 = exp[1][0]
189 :     e_20 = exp[2][0]
190 :     e_01 = exp[0][1]
191 :     e_11 = exp[1][1]
192 :     e_21 = exp[2][1]
193 :     e_02 = exp[0][2]
194 :     e_12 = exp[1][2]
195 :     e_22 = exp[2][2]
196 :     def diff_d3(e):
197 :     return [diff(e,x), diff(e,y), diff(e,z)]
198 :     f_00 = diff_d3(e_00)
199 :     f_01 = diff_d3(e_01)
200 :     f_02 = diff_d3(e_02)
201 :     f_10 = diff_d3(e_10)
202 :     f_11 = diff_d3(e_11)
203 :     f_12 = diff_d3(e_12)
204 :     f_20 = diff_d3(e_20)
205 :     f_21 = diff_d3(e_21)
206 :     f_22 = diff_d3(e_22)
207 :     return [[f_00,f_01,f_02],[f_10,f_11,f_12],[f_20,f_21,f_22]]
208 :     else:
209 :     raise (" type not supported for jacob")
210 : cchiw 3939 else:
211 : cchiw 3998 raise "unsupported type for jacob"
212 : cchiw 3939
213 : cchiw 3998 #evaluate norm
214 :     def fn_norm(fld, dim):
215 :     exp = field.get_data(fld)
216 :     ityp = field.get_ty(fld)
217 :     dim = field.get_dim(fld)
218 :     #print " exp1: ",exp1," exp2: ",exp2
219 :     # vectors
220 :     def iter (es):
221 :     sum = 0
222 :     for i in es:
223 : cchiw 4334 t=i*i
224 :     #print "t",t
225 :     sum+=t
226 : cchiw 4210 #print "\nsum",sum
227 : cchiw 4250 rtn = sqrt(sum)
228 : cchiw 4210 #print "\nrtn",rtn
229 : cchiw 4188 return rtn
230 : cchiw 3998 if(field.is_Scalar(fld)):
231 :     [] = fty.get_shape(ityp)
232 : cchiw 4334 #print "scalar exp:",exp
233 : cchiw 4343 t =(sqrt(exp*exp))
234 : cchiw 4334 #print "t",t
235 :     return t
236 : cchiw 3998 elif(field.is_Vector(fld)):
237 :     [n] = fty.get_shape(ityp)
238 :     rtn = []
239 :     for i in range(n):
240 :     rtn.append(exp[i])
241 :     return iter(rtn)
242 :     elif(field.is_Matrix(fld)):
243 :     [n, m] = fty.get_shape(ityp)
244 :     rtn = []
245 :     for i in range(n):
246 :     for j in range(m):
247 :     rtn.append(exp[i][j])
248 :     return iter(rtn)
249 :     elif(field.is_Ten3(fld)):
250 :     [n, m, p] = fty.get_shape(ityp)
251 :     rtn = []
252 :     for i in range(n):
253 :     for j in range(m):
254 :     for k in range(p):
255 :     rtn.append(exp[i][j][k])
256 :     return iter(rtn)
257 :     else:
258 :     raise "unsupported type for norm"
259 : cchiw 3939
260 : cchiw 3998 #evaluate norm
261 :     def fn_normalize(fld, dim):
262 :     exp = field.get_data(fld)
263 :     ityp = field.get_ty(fld)
264 :     dim = field.get_dim(fld)
265 :     #print " exp1: ",exp1," exp2: ",exp2
266 :     norm = fn_norm(fld, dim)
267 :     if(field.is_Scalar(fld)):
268 :     #print "scal",exp
269 :     return exp
270 :     elif(field.is_Vector(fld)):
271 :     [n] = fty.get_shape(ityp)
272 :     rtn = []
273 :     for i in range(n):
274 :     rtn.append(exp[i]/norm)
275 :     #print "vec",rtn
276 :     return rtn
277 :     elif(field.is_Matrix(fld)):
278 :     [n, m] = fty.get_shape(ityp)
279 :     rtn = []
280 :     for i in range(n):
281 :     rtni = []
282 :     for j in range(m):
283 :     rtni.append(exp[i][j]/norm)
284 :     rtn.append(rtni)
285 :     #print "matrix:",rtn
286 :     return rtn
287 :     elif(field.is_Ten3(fld)):
288 :     [n, m, p] = fty.get_shape(ityp)
289 :     rtn = []
290 :     for i in range(n):
291 :     rtni = []
292 :     for j in range(m):
293 :     rtnj = []
294 :     for k in range(p):
295 :     rtnj.append(exp[i][j][k]/norm)
296 :     rtni.append( rtnj)
297 :     rtn.append( rtni)
298 :     #print "ten3",rtn
299 :     return rtn
300 :     else:
301 :     raise "unsupported type for norm"
302 :    
303 : cchiw 4237 #[0]
304 :     def fn_slicev0(fld1):
305 :     exp1 = field.get_data(fld1)
306 :     ityp1 = field.get_ty(fld1)
307 :     return exp1[0]
308 :    
309 :    
310 : cchiw 3998 #evaluate slice
311 : cchiw 4237 #[1]
312 :     def fn_slicev1(fld1):
313 : cchiw 3998 exp1 = field.get_data(fld1)
314 :     ityp1 = field.get_ty(fld1)
315 : cchiw 4237 return exp1[1]
316 :    
317 :     #evaluate slice
318 :     #[1,:]
319 :     def fn_slicem0(fld1):
320 :     exp1 = field.get_data(fld1)
321 :     ityp1 = field.get_ty(fld1)
322 : cchiw 3998 rtn=[]
323 :     if(fty.is_Matrix(ityp1)):
324 :     [n2,n3] = fty.get_shape(ityp1)
325 : cchiw 4237 for i in range(n3):
326 :     rtn.append(exp1[1][i])
327 : cchiw 3998 return rtn
328 :     else:
329 :     raise "unsupported type for slice"
330 :    
331 : cchiw 4237 #evaluate slice
332 :     #[:,0]
333 :     def fn_slicem1(fld1):
334 :     exp1 = field.get_data(fld1)
335 :     ityp1 = field.get_ty(fld1)
336 :     rtn=[]
337 :     if(fty.is_Matrix(ityp1)):
338 :     [n2,n3] = fty.get_shape(ityp1)
339 :     for i in range(n2):
340 :     rtn.append(exp1[i][0])
341 :     return rtn
342 :     else:
343 :     raise "unsupported type for slice"
344 :    
345 : cchiw 4243 #evaluate slice
346 :     #[:,1,:]
347 :     def fn_slicet0(fld1):
348 :     exp1 = field.get_data(fld1)
349 :     ityp1 = field.get_ty(fld1)
350 :     rtn=[]
351 :     if(fty.is_Ten3(ityp1)):
352 :     [n1, n2,n3] = fty.get_shape(ityp1)
353 :     for i in range(n1):
354 :     rtnj=[]
355 :     for j in range(n3):
356 :     rtnj.append(exp1[i][1][j])
357 :     rtn.append(rtnj)
358 :     return rtn
359 :     else:
360 :     raise "unsupported type for slice"
361 : cchiw 4237
362 : cchiw 4243 #evaluate slice
363 :     #[1,0,:]
364 :     def fn_slicet1(fld1):
365 :     exp1 = field.get_data(fld1)
366 :     ityp1 = field.get_ty(fld1)
367 :     rtn=[]
368 :     if(fty.is_Ten3(ityp1)):
369 :     [n1, n2, n3] = fty.get_shape(ityp1)
370 :     for i in range(n3):
371 :     rtn.append(exp1[1][0][i])
372 :     return rtn
373 :     else:
374 :     raise "unsupported type for slice"
375 :    
376 :    
377 :    
378 :    
379 : cchiw 3998 #evaluate trace
380 :     def fn_trace(fld):
381 :     exp = field.get_data(fld)
382 :     ityp = field.get_ty(fld)
383 :     rtn=[]
384 :     if(field.is_Matrix(fld)):
385 :     [n, m] = fty.get_shape(ityp)
386 :     if (n!=m):
387 :     raise Exception("matrix is not identitical")
388 :     rtn = exp[0][0]+exp[1][1]
389 :     if(n==2):
390 :     return rtn
391 :     elif(n==3):
392 :     return rtn+exp[2][2]
393 :     else:
394 :     raise "unsupported matrix field"
395 :     else:
396 :     raise "unsupported trace"
397 :    
398 :     #evaluate transpose
399 :     def fn_transpose(fld):
400 :     exp = field.get_data(fld)
401 :     ityp = field.get_ty(fld)
402 :     if(field.is_Matrix(fld)):
403 :     [n, m] = fty.get_shape(ityp)
404 :     rtn = []
405 :     for i in range(n):
406 :     rtni = []
407 :     for j in range(m):
408 :     rtni.append(exp[j][i])
409 :     rtn.append(rtni)
410 :     return rtn
411 :     else:
412 :     raise "unsupported transpose"
413 :    
414 : cchiw 4247 def fn_doubledot(fld1, fld2):
415 :     exp1 = field.get_data(fld1)
416 :     exp2 = field.get_data(fld2)
417 :     ityp1 = field.get_ty(fld1)
418 : cchiw 4308 if(field.is_Matrix(fld1)):
419 :     rtn = 0
420 :     [n, m] = fty.get_shape(ityp1)
421 :     for i in range(n):
422 : cchiw 4247 for j in range(m):
423 : cchiw 4248 rtn+=exp1[i][j]*exp2[i][j]
424 : cchiw 4308 return rtn
425 :     elif(field.is_Ten3(fld1)):
426 :     rtn = []
427 :     [n, m, l] = fty.get_shape(ityp1)
428 :     for i in range(n):
429 :     for j in range(m):
430 :     for k in range(l):
431 :     rtn.append(exp1[i][j][k]*exp2[i][j][k])
432 :     return rtn
433 :     else:
434 :     raise "unsupported double dot"
435 : cchiw 4247
436 :    
437 :    
438 : cchiw 3998 #evaluate det
439 :     def fn_det(fld):
440 :     exp = field.get_data(fld)
441 :     ityp = field.get_ty(fld)
442 :     rtn=[]
443 :     if(field.is_Matrix(fld)):
444 :     [n, m] = fty.get_shape(ityp)
445 :     if (n!=m):
446 :     raise Exception("matrix is not identitical")
447 :     a = exp[0][0]
448 :     d = exp[1][1]
449 :     c = exp[1][0]
450 :     b = exp[0][1]
451 :     if(n==2):
452 :     return a*d-b*c
453 :     elif(n==3):
454 :     a = exp[0][0]
455 :     b = exp[0][1]
456 :     c = exp[0][2]
457 :     d = exp[1][0]
458 :     e = exp[1][1]
459 :     f = exp[1][2]
460 :     g = exp[2][0]
461 :     h = exp[2][1]
462 :     i = exp[2][2]
463 :     return a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)
464 :     else:
465 :     raise "unsupported matrix field"
466 :     else:
467 :     raise "unsupported trace"
468 :    
469 : cchiw 4411 #evaluate det
470 :     def fn_inverse(fld):
471 :     exp = field.get_data(fld)
472 :     ityp = field.get_ty(fld)
473 :     rtn=[]
474 :     if(field.is_Matrix(fld)):
475 :     [n, m] = fty.get_shape(ityp)
476 :     if (n!=m):
477 :     raise Exception("matrix is not identitical")
478 :     a = exp[0][0]
479 :     d = exp[1][1]
480 :     c = exp[1][0]
481 :     b = exp[0][1]
482 :     if(n==2):
483 :     denom = a*d-b*c
484 :     return [[d/denom,-b/denom],[-c/denom, a/denom]]
485 :     elif(n==3):
486 :     a = exp[0][0]
487 :     b = exp[0][1]
488 :     c = exp[0][2]
489 :     d = exp[1][0]
490 :     e = exp[1][1]
491 :     f = exp[1][2]
492 :     g = exp[2][0]
493 :     h = exp[2][1]
494 :     i = exp[2][2]
495 :     denom = a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)
496 :     num1 = [(e*i - f*h), -(b*i - c*h), (b*f - c*e)]
497 :     num2 = [-(d*i - f*g), (a*i -c*g), -(a*f - c*d)]
498 :     num3 = [(d*h - e*g), -(a*h - b*g),(a*e - b*d)]
499 :     num =[num1,num2,num3]
500 :    
501 :     rtn = []
502 :     for i in range(n):
503 :     tmp = []
504 :     for j in range(n):
505 :     tmp.append(num[i][j]/denom)
506 :     rtn.append(tmp)
507 :     return rtn
508 :     else:
509 :     raise "unsupported matrix field"
510 :     else:
511 :     raise "unsupported trace"
512 : cchiw 3998
513 : cchiw 4411
514 :    
515 :    
516 :    
517 : cchiw 3939 #evaluate outer product
518 :     def fn_outer(fld1, fld2):
519 :     exp1 = field.get_data(fld1)
520 :     ityp1 = field.get_ty(fld1)
521 :     exp2 = field.get_data(fld2)
522 :     ityp2 = field.get_ty(fld2)
523 : cchiw 4268 ashape = fty.get_shape(ityp1)
524 :     bshape = fty.get_shape(ityp2)
525 :     x= "ashape", ashape, "bshape", bshape
526 : cchiw 3946 #print "exp1",exp1,"ityp1",ityp1.name,"-length",len(exp1)
527 :     #print "exp2",exp2,"ityp2",ityp2.name,"-length",len(exp2)
528 : cchiw 4268 rtn = []
529 : cchiw 3939 if(fty.is_Vector(ityp1)):
530 : cchiw 4268 [n1] = fty.get_shape(ityp1)
531 : cchiw 3939 if(fty.is_Vector(ityp2)):
532 :     #both vectors
533 : cchiw 4268 [n2] = fty.get_shape(ityp2)
534 : cchiw 4456 x= "\n outer made shape:"+(str(n1)+","+str(n2))
535 :     write_outer(x)
536 : cchiw 3939 for i in range(n1):
537 : cchiw 3946 tmpI = []
538 : cchiw 3939 for j in range(n2):
539 : cchiw 3946 tmpI.append(exp1[i]*exp2[j])
540 :     rtn.append(tmpI)
541 : cchiw 3939 return rtn
542 :     elif(fty.is_Matrix(ityp2)):
543 :     [n2,n3] = fty.get_shape(ityp2)
544 : cchiw 4456 x= "\n outer made shape:"+(str(n1)+","+str(n2)+","+str(n3))
545 :     write_outer(x)
546 : cchiw 3939 for i in range(n1):
547 : cchiw 3946 tmpI = []
548 : cchiw 3939 for j in range(n2):
549 : cchiw 3946 tmpJ = []
550 : cchiw 3939 for k in range(n3):
551 : cchiw 3946 tmpJ.append(exp1[i]*exp2[j][k])
552 :     tmpI.append(tmpJ)
553 :     rtn.append(tmpI)
554 : cchiw 3939 return rtn
555 :     elif(fty.is_Matrix(ityp1)):
556 :     [n1,n2] = fty.get_shape(ityp1)
557 :     if(fty.is_Vector(ityp2)):
558 : cchiw 4268 [n3] = fty.get_shape(ityp2)
559 : cchiw 4456 x= "\n outer made shape:"+(str(n1)+","+str(n2)+","+str(n3))
560 :     write_outer(x)
561 : cchiw 3939 for i in range(n1):
562 : cchiw 3946 tmpI = []
563 : cchiw 3939 for j in range(n2):
564 : cchiw 3946 tmpJ = []
565 : cchiw 3939 for k in range(n3):
566 : cchiw 3946 tmpJ.append(exp1[i][j]*exp2[k])
567 :     tmpI.append(tmpJ)
568 :     rtn.append(tmpI)
569 : cchiw 3939 return rtn
570 :     elif(fty.is_Matrix(ityp2)):
571 : cchiw 4268 [n3, n4] = fty.get_shape(ityp2)
572 : cchiw 4456 x= "\n outer made shape:"+(str(n1)+","+str(n4))
573 :     write_outer(x)
574 : cchiw 3939 for i in range(n1):
575 : cchiw 3946 tmpI = []
576 : cchiw 3939 for j in range(n2):
577 : cchiw 3946 tmpJ = []
578 : cchiw 3939 for k in range(n3):
579 : cchiw 3946 tmpK = []
580 : cchiw 3939 for l in range(n4):
581 : cchiw 3946 tmpK.append(exp1[i][j]*exp2[k][l])
582 :     tmpJ.append(tmpK)
583 :     tmpI.append(tmpJ)
584 :     rtn.append(tmpI)
585 : cchiw 3939 return rtn
586 :     else:
587 :     raise "outer product is not supported"
588 :     else:
589 :     raise "outer product is not supported"
590 :     #evaluate inner product
591 :     def fn_inner(fld1, fld2):
592 :     exp1 = field.get_data(fld1)
593 :     ityp1 = field.get_ty(fld1)
594 :     exp2 = field.get_data(fld2)
595 :     ityp2 = field.get_ty(fld2)
596 : cchiw 4268 ashape = fty.get_shape(ityp1)
597 :     bshape = fty.get_shape(ityp2)
598 :     x= "ashape", ashape, "bshape", bshape
599 : cchiw 3939 if(fty.is_Vector(ityp1)):
600 : cchiw 4268 [a1] = fty.get_shape(ityp1)
601 : cchiw 3939 if(fty.is_Vector(ityp2)):
602 :     #length of vetors
603 :     rtn=0
604 : cchiw 4269 [b1] = fty.get_shape(ityp2)
605 : cchiw 4268 if(a1!=b1):
606 :     raise x
607 :     for s in range(a1):
608 :     rtn += exp1[s]*exp2[s]
609 : cchiw 3939 return rtn
610 :     elif(fty.is_Matrix(ityp2)):
611 : cchiw 4268 [b1,b2] = fty.get_shape(ityp2)
612 : cchiw 3939 rtn=[]
613 : cchiw 4268 if(a1!=b1):
614 :     raise x
615 :     for i in range(b2):
616 : cchiw 3939 sumrtn=0
617 : cchiw 4268 for s in range(a1):
618 :     sumrtn += exp1[s]*exp2[s][i]
619 : cchiw 3939 rtn.append(sumrtn)
620 :     return rtn
621 :     elif(fty.is_Ten3(ityp2)):
622 : cchiw 4268 [b1, b2, b3] = fty.get_shape(ityp2)
623 : cchiw 3946 rtn = []
624 : cchiw 4268 if(a1!=b1):
625 :     raise x
626 :     for i in range(b2):
627 : cchiw 3946 tmpJ = []
628 : cchiw 4268 for j in range(b3):
629 : cchiw 3939 sumrtn=0
630 : cchiw 4268 for s in range(a1):
631 :     sumrtn += exp1[s]*exp2[s][i][j]
632 : cchiw 3946 tmpJ.append(sumrtn)
633 :     rtn.append(tmpJ)
634 :     return rtn
635 : cchiw 3939 else:
636 :     raise "inner product is not supported"
637 :     elif(fty.is_Matrix(ityp1)):
638 : cchiw 4268 [a1,a2] = fty.get_shape(ityp1)
639 : cchiw 3939 if(fty.is_Vector(ityp2)):
640 : cchiw 4268 [b1] = fty.get_shape(ityp2)
641 :     if(a2!=b1):
642 :     raise x
643 : cchiw 3939 rtn=[]
644 : cchiw 4268 for i in range(a1):
645 : cchiw 3939 sumrtn=0
646 : cchiw 4268 for s in range(a2):
647 :     sumrtn += exp1[i][s]*exp2[s]
648 : cchiw 3939 rtn.append(sumrtn)
649 :     return rtn
650 : cchiw 4308 elif(fty.is_Matrix(ityp2)):
651 : cchiw 4261 [b1,b2] = fty.get_shape(ityp2)
652 :     rtn=[]
653 : cchiw 4268 if(a2!=b1):
654 :     raise x
655 : cchiw 4261 for i in range(a1):
656 :     rtnj = []
657 :     for j in range(b2):
658 :     sumrtn=0
659 :     for s in range(a2):
660 : cchiw 4268 sumrtn += exp1[i][s]*exp2[s][j]
661 : cchiw 4261 rtnj.append(sumrtn)
662 :     rtn.append(rtnj)
663 :     return rtn
664 : cchiw 4308 elif(fty.is_Ten3(ityp2)):
665 :     [b1,b2, b3] = fty.get_shape(ityp2)
666 :     rtn=[]
667 :     if(a2!=b1):
668 :     raise x
669 :     for i in range(a1):
670 :     rtnj = []
671 :     for j in range(b2):
672 :     rtnk = []
673 :     for k in range(b3):
674 :     sumrtn=0
675 :     for s in range(a2):
676 :     sumrtn += exp1[i][s]*exp2[s][j][k]
677 :     rtnk.append(sumrtn)
678 :     rtnj.append(rtnk)
679 :     rtn.append(rtnj)
680 :     return rtn
681 :    
682 : cchiw 3939 else:
683 :     raise "inner product is not supported"
684 :     elif(fty.is_Ten3(ityp1)):
685 : cchiw 4268 [a1,a2, a3] = ashape
686 : cchiw 3939 if(fty.is_Vector(ityp2)):
687 : cchiw 4268 [b1] = bshape
688 :     if(a3!=b1):
689 :     raise x
690 : cchiw 3939 rtn=[]
691 : cchiw 4268 for i in range(a1):
692 : cchiw 3946 tmpI=[]
693 : cchiw 4268 for j in range(a2):
694 : cchiw 3939 sumrtn=0
695 : cchiw 4268 for s in range(a3):
696 :     sumrtn += exp1[i][j][s]*exp2[s]
697 : cchiw 3946 tmpI.append(sumrtn)
698 :     rtn.append(tmpI)
699 : cchiw 3939 return rtn
700 : cchiw 4268 if(fty.is_Matrix(ityp2)):
701 :     [b1,b2] = bshape
702 :     if(a3!=b1):
703 :     raise x
704 :     rtn=[]
705 :     for i in range(a1):
706 :     tmpI=[]
707 :     for j in range(a2):
708 :     tmpJ = []
709 :     for k in range(b2):
710 :     sumrtn=0
711 :     for s in range(a3):
712 :     sumrtn += exp1[i][j][s]*exp2[s][k]
713 :     tmpJ.append(sumrtn)
714 :     tmpI.append(tmpJ)
715 :     rtn.append(tmpI)
716 :     return rtn
717 : cchiw 3939 else:
718 :     raise "inner product is not supported"
719 :     else:
720 :     raise "inner product is not supported"
721 :    
722 :    
723 : cchiw 3915 # *************************** generic apply operators ***************************
724 :     #unary operator on a vector
725 : cchiw 3939 def applyToVector(vec, unary):
726 : cchiw 3915 rtn = []
727 :     for v in vec:
728 :     rtn.append(unary(v))
729 :     return rtn
730 :     #binary operator on a vector
731 : cchiw 3939 def applyToVectors(vecA, vecB, binary):
732 : cchiw 3915 rtn = []
733 :     for (a,b) in zip(vecA,vecB):
734 :     x= binary(a,b)
735 :     rtn.append(x)
736 :     return rtn
737 : cchiw 3874
738 : cchiw 3946 def applyToM(vec, unary):
739 :     rtn = []
740 :     for i in vec:
741 :     tmpI = []
742 :     for v in i:
743 :     tmpI.append(unary(v))
744 :     rtn.append(tmpI)
745 :     return rtn
746 :    
747 : cchiw 4230 def applyToMs(vecA,vecB, unary):
748 :     rtn = []
749 :     for (a,b) in zip(vecA,vecB):
750 :     tmpI = []
751 :     for (u,v) in zip(a,b):
752 :     tmpI.append(unary(u, v))
753 :     rtn.append(tmpI)
754 :     return rtn
755 :    
756 :    
757 : cchiw 3946 def applyToT3(vec, unary):
758 :     rtn = []
759 :     for i in vec:
760 :     tmpI = []
761 :     for j in i:
762 :     tmpJ = []
763 :     for v in j:
764 :     tmpJ.append(unary(v))
765 :     tmpI.append(tmpJ)
766 :     rtn.append(tmpI)
767 :     return rtn
768 :    
769 : cchiw 4308 def applyToT3s(vecA, vecB, unary):
770 :     rtn = []
771 :     for (a,b) in zip(vecA, vecB):
772 :     tmpI = []
773 :     for (i,j) in zip(a,b):
774 :     tmpJ = []
775 :     for (u,v) in zip(i,j):
776 :     tmpJ.append(unary(u, v))
777 :     tmpI.append(tmpJ)
778 :     rtn.append(tmpI)
779 :     return rtn
780 :    
781 :    
782 : cchiw 3915 # *************************** apply to scalars or vectors ***************************
783 :     #apply operators to expression
784 :     # return output types and expression
785 :     # unary operator
786 :     # exp: scalar types
787 : cchiw 3874
788 : cchiw 3939 def applyToExp_U_S(fn_name, fld):
789 :     exp = field.get_data(fld)
790 :     dim = field.get_dim(fld)
791 :     #print fn_name
792 : cchiw 4338 if(op_copy==fn_name): #probing
793 : cchiw 3939 return exp
794 : cchiw 3915 elif(op_negation==fn_name): #negation
795 : cchiw 3939 return fn_negation(exp)
796 : cchiw 4250 elif(op_cosine==fn_name): #cosine
797 :     return cos(exp)
798 :     elif(op_sine==fn_name): #sine
799 :     return sin(exp)
800 :     elif(op_tangent==fn_name): # tangent
801 :     return tan(exp)
802 :     elif(op_atangent==fn_name): #atan
803 :     return atan(exp)
804 : cchiw 4252 elif(op_gradient==fn_name):
805 :     return fn_grad(exp, dim)
806 : cchiw 4444 elif(op_asine==fn_name): #asine mag(x)<=1
807 : cchiw 4455 frac = 0.01*exp
808 : cchiw 4444 return asin(frac)
809 :     elif(op_acosine==fn_name): #acos mag(x)<=1
810 : cchiw 4455 frac = 0.01*exp
811 : cchiw 4444 return acos(frac)
812 :     elif(op_sqrt==fn_name): #sqrt
813 :     # gets norm first to make sure value is positive.
814 :     norm = sqrt(exp*exp)
815 :     return sqrt(norm)
816 : cchiw 3915 else:
817 : cchiw 3939 raise Exception("unsupported unary operator on scalar field:"+ fn_name.name)
818 : cchiw 3874
819 : cchiw 3939 # unary operator
820 :     # exp: vector types
821 :     def applyToExp_U_V(fn_name, fld):
822 :     exp = field.get_data(fld)
823 : cchiw 4338 if(op_copy==fn_name): #probing
824 : cchiw 3939 return exp
825 :     elif(op_negation==fn_name): #negation
826 :     return applyToVector(exp, fn_negation)
827 :     elif(op_divergence==fn_name):
828 :     return fn_divergence(fld)
829 :     elif(op_curl==fn_name):
830 :     return fn_curl(fld)
831 :     elif(op_jacob==fn_name): #jacob
832 :     return fn_jacob(fld)
833 : cchiw 4237 elif(op_slicev0==fn_name) :
834 :     return fn_slicev0(fld)
835 :     elif(op_slicev1==fn_name):
836 :     return fn_slicev1(fld)
837 : cchiw 4252 elif(op_grad==op_grad):
838 :     return fn_grad(fld)
839 : cchiw 3939 else:
840 :     raise Exception("unsupported unary operator:"+ fn_name.name)
841 : cchiw 3874
842 : cchiw 3946 def applyToExp_U_M(fn_name, fld):
843 :     exp = field.get_data(fld)
844 : cchiw 4338 if(op_copy==fn_name): #probing
845 : cchiw 3946 return exp
846 :     elif(op_negation==fn_name): #negation
847 :     return applyToM(exp, fn_negation)
848 :     elif(op_jacob==fn_name): #jacob
849 : cchiw 4411 #print "app u-m"
850 :     x = fn_jacob(fld)
851 :     #print "x", x
852 :     return x
853 : cchiw 4237 elif(op_slicem0==fn_name) :
854 :     return fn_slicem0(fld)
855 :     elif(op_slicem1==fn_name):
856 :     return fn_slicem1(fld)
857 : cchiw 3998 elif(op_trace == fn_name):
858 :     return fn_trace(fld)
859 :     elif(op_transpose==fn_name):
860 :     return fn_transpose(fld)
861 :     elif(op_det==fn_name):
862 :     return fn_det(fld)
863 : cchiw 4411 elif(op_inverse == fn_name):
864 :     return fn_inverse(fld)
865 : cchiw 3946 else:
866 :     raise Exception("unsupported unary operator:"+ fn_name.name)
867 :    
868 : cchiw 3939 def applyToExp_U_T3(fn_name, fld):
869 :     exp = field.get_data(fld)
870 : cchiw 4338 if(op_copy==fn_name): #probing
871 : cchiw 3939 return exp
872 :     elif(op_negation==fn_name): #negation
873 : cchiw 3946 return applyToT3(exp, fn_negation)
874 : cchiw 3939 elif(op_jacob==fn_name): #jacob
875 :     return fn_jacob(fld)
876 : cchiw 4243 elif(op_slicet0==fn_name) :
877 :     return fn_slicet0(fld)
878 :     elif(op_slicet1==fn_name):
879 :     return fn_slicet1(fld)
880 : cchiw 3939 else:
881 :     raise Exception("unsupported unary operator:"+ fn_name.name)
882 :    
883 : cchiw 4308 # binary, args do not need to have the same shape
884 :     def applyToExp_B_rest(e):
885 :     fn_name=e.opr
886 :     (fld1,fld2) = apply.get_binary(e)
887 :     exp1 = field.get_data(fld1)
888 :     exp2 = field.get_data(fld2)
889 :     if(op_outer==fn_name):
890 :     return fn_outer(fld1, fld2)
891 :     elif(op_inner==fn_name):
892 :     return fn_inner(fld1, fld2)
893 :     elif(op_scale==fn_name): #scaling
894 :     return fn_scaling(fld1,fld2)
895 :     else:
896 :     raise Exception("unsupported unary operator:"+fn_name.name)
897 :    
898 : cchiw 3915 # binary operator
899 :     # exp: scalar types
900 :     def applyToExp_B_S(e):
901 :     fn_name=e.opr
902 : cchiw 4308
903 : cchiw 3939 (fld1,fld2) = apply.get_binary(e)
904 :     exp1 = field.get_data(fld1)
905 :     exp2 = field.get_data(fld2)
906 :     #print fn_name
907 : cchiw 3915 if(op_add==fn_name):#addition
908 : cchiw 3939 return fn_add(exp1,exp2)
909 : cchiw 3915 elif(op_subtract==fn_name):#subtract
910 : cchiw 3939 return fn_subtract(exp1,exp2)
911 : cchiw 4243 elif(op_modulate==fn_name):#modulate
912 :     return fn_modulate(exp1,exp2)
913 : cchiw 3939 elif(op_scale==fn_name): #scaling
914 :     return fn_scaling(fld1,fld2)
915 : cchiw 3998 elif(op_division==fn_name): #division
916 :     return fn_division(fld1,fld2)
917 : cchiw 3915 else:
918 : cchiw 3939 raise Exception("unsupported binary operator on scalar fields:"+ fn_name.name)
919 : cchiw 3874
920 : cchiw 3865
921 : cchiw 3915 # binary operator
922 : cchiw 3939 # args have the same shape
923 : cchiw 3915 def applyToExp_B_V(e):
924 : cchiw 4308
925 : cchiw 3915 fn_name=e.opr
926 : cchiw 3939 (fld1,fld2) = apply.get_binary(e)
927 :     exp1 = field.get_data(fld1)
928 :     exp2 = field.get_data(fld2)
929 : cchiw 3915 if(op_add==fn_name):#addition
930 : cchiw 3939 return applyToVectors(exp1, exp2, fn_add)
931 : cchiw 3915 elif(op_subtract==fn_name):#subtract
932 : cchiw 3939 return applyToVectors(exp1, exp2, fn_subtract)
933 : cchiw 4243 elif(op_modulate==fn_name):#modulate
934 :     return applyToVectors(exp1,exp2 ,fn_modulate)
935 : cchiw 3939 elif(op_cross==fn_name):
936 :     return fn_cross(fld1, fld2)
937 : cchiw 3915 else:
938 : cchiw 4308 return applyToExp_B_rest(e)
939 : cchiw 3915
940 : cchiw 4230 def applyToExp_B_M(e):
941 : cchiw 4261 print "b. m"
942 : cchiw 4230 fn_name=e.opr
943 :     (fld1,fld2) = apply.get_binary(e)
944 :     exp1 = field.get_data(fld1)
945 :     exp2 = field.get_data(fld2)
946 :     if(op_add==fn_name):#addition
947 :     return applyToMs(exp1, exp2, fn_add)
948 :     elif(op_subtract==fn_name):#subtract
949 :     return applyToMs(exp1, exp2, fn_subtract)
950 : cchiw 4243 elif(op_modulate==fn_name):#modulate
951 :     return applyToMs(exp1,exp2,fn_modulate)
952 : cchiw 4230 else:
953 : cchiw 4308 return applyToExp_B_rest(e)
954 : cchiw 3939
955 : cchiw 4308 def applyToExp_B_T3(e):
956 :     print "b. m"
957 :     fn_name=e.opr
958 :     (fld1,fld2) = apply.get_binary(e)
959 :     exp1 = field.get_data(fld1)
960 :     exp2 = field.get_data(fld2)
961 :     if(op_add==fn_name):#addition
962 :     return applyToT3s(exp1, exp2, fn_add)
963 :     elif(op_subtract==fn_name):#subtract
964 :     return applyToT3s(exp1, exp2, fn_subtract)
965 :     elif(op_modulate==fn_name):#modulate
966 :     return applyToT3s(exp1,exp2,fn_modulate)
967 :     else:
968 :     return applyToExp_B_rest(e)
969 : cchiw 4230
970 : cchiw 4308
971 : cchiw 3915 # *************************** unary/binary operators ***************************
972 :     def unary(e):
973 : cchiw 3939 #apply.toStr(e)
974 :     fld =apply.get_unary(e)
975 :     fn_name=e.opr
976 :     exp = field.get_data(fld)
977 :     dim = field.get_dim(fld)
978 : cchiw 3998 if(op_norm==fn_name):#norm
979 :     return fn_norm(fld, dim)
980 : cchiw 4308 elif(op_normalize==fn_name):#normalize
981 : cchiw 3998 x= fn_normalize(fld, dim)
982 :     return x
983 :     elif (field.is_Scalar(fld)): # input is a scalar field
984 : cchiw 3939 return applyToExp_U_S(fn_name, fld)
985 :     elif(field.is_Vector(fld)): # input is a vector field
986 :     return applyToExp_U_V(fn_name, fld)
987 : cchiw 3946 elif(field.is_Matrix(fld)): # input is a vector field
988 :     return applyToExp_U_M(fn_name, fld)
989 : cchiw 3915 else:
990 : cchiw 3939 return applyToExp_U_T3(fn_name, fld)
991 : cchiw 3915
992 :     def binary(e):
993 : cchiw 3939 (f, g) =apply.get_binary(e)
994 : cchiw 3998 fn_name = e.opr
995 :     # type is checked elsewhere or does not matter
996 :     if(op_division==fn_name): #division
997 :     return fn_division(f, g)
998 : cchiw 4308 elif(op_doubledot==fn_name):#op_doubledot
999 :     return fn_doubledot (f, g)
1000 :     elif(op_outer==fn_name):
1001 :     return fn_outer(f, g)
1002 :     elif(op_inner==fn_name):
1003 :     return fn_inner(f, g)
1004 :     elif(op_scale==fn_name): #scaling
1005 :     return fn_scaling(f,g)
1006 : cchiw 3998 elif (field.is_Scalar(f) and field.is_Scalar(g)): # input is a scalar field
1007 : cchiw 3915 return applyToExp_B_S(e)
1008 : cchiw 4308 elif (field.is_Vector(f) and field.is_Vector(g)):
1009 : cchiw 3939 return applyToExp_B_V(e)
1010 : cchiw 4308 elif (field.is_Matrix(f) and field.is_Matrix(g)):
1011 :     return applyToExp_B_M(e)
1012 :     elif (field.is_Ten3(f) and field.is_Ten3(g)):
1013 :     return applyToExp_B_T3(e)
1014 : cchiw 3915 else:
1015 : cchiw 4308 return applyToExp_B_rest(e)
1016 : cchiw 3915
1017 : cchiw 3939 def applyUnaryOnce(oexp_inner,app_inner,app_outer):
1018 : cchiw 3998 #print "applyUnaryOnce"
1019 : cchiw 3939 #apply.toStr(app_inner)
1020 :     oty_inner = apply.get_oty(app_inner)
1021 :     oty_outer = apply.get_oty(app_outer)
1022 :     opr_outer = app_outer.opr
1023 : cchiw 3946 #print "oexp_inner",oexp_inner,"opr_outer",opr_outer.name
1024 : cchiw 3939 lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
1025 : cchiw 3946 app_tmp = apply("tmp", opr_outer, lhs_tmp, None, oty_outer, true, true)
1026 : cchiw 3939 oexp_tmp =unary(app_tmp)
1027 : cchiw 3946 #print " oexp_tmp", oexp_tmp
1028 : cchiw 3939 return (oty_outer, oexp_tmp)
1029 :    
1030 : cchiw 3946 def applyBinaryOnce(oexp_inner,app_inner,app_outer,rhs):
1031 :     oty_inner = apply.get_oty(app_inner)
1032 :     oty_outer = apply.get_oty(app_outer)
1033 :     opr_outer = app_outer.opr
1034 :    
1035 :     lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
1036 :    
1037 :     app_tmp = apply("tmp", opr_outer, lhs_tmp, rhs, oty_outer, true,true)
1038 :     oexp_tmp =binary(app_tmp)
1039 :     return (oty_outer, oexp_tmp)
1040 : cchiw 3939
1041 : cchiw 3946
1042 : cchiw 3915 # operators with scalar field and vector field
1043 :     def sort(e):
1044 : cchiw 3939 #apply.toStr(e)
1045 :     arity = apply.get_arity(e)
1046 : cchiw 3946 if(e.isrootlhs): # is root
1047 : cchiw 3939 #print "sort is a root"
1048 :     oty = apply.get_oty(e)
1049 :     if (arity ==1):
1050 :     return (oty, unary(e))
1051 :     elif (arity ==2): # binary operator
1052 :     return (oty, binary(e))
1053 :     else:
1054 :     raise Exception ("arity is not supported: "+str(arity_outer))
1055 : cchiw 3915 else:
1056 : cchiw 3939 app_outer = e
1057 :     arity_outer = arity
1058 : cchiw 3946 #print "app_outer",app_outer.opr.name
1059 : cchiw 3939 if (arity_outer==1): #assuming both arity
1060 :     app_inner = apply.get_unary(app_outer)
1061 : cchiw 3946 #print "outer(1) app_inner:",app_inner.opr.name
1062 : cchiw 3939 arity_inner= app_inner.opr.arity
1063 :     if (arity_inner==1):
1064 :     oexp_inner = unary(app_inner)
1065 : cchiw 3946
1066 : cchiw 3939 (oty_outer, oexp_tmp) = applyUnaryOnce(oexp_inner ,app_inner, app_outer)
1067 : cchiw 3946
1068 : cchiw 3939 return (oty_outer, oexp_tmp)
1069 :     elif(arity_inner==2):
1070 :     oexp_inner = binary(app_inner)
1071 :     (oty_outer, oexp_tmp) = applyUnaryOnce(oexp_inner ,app_inner, app_outer)
1072 : cchiw 3946
1073 : cchiw 3939 return (oty_outer, oexp_tmp)
1074 :     else:
1075 :     raise Exception ("arity is not supported: "+str(arity_outer))
1076 : cchiw 3946 elif (arity_outer==2):
1077 :     (app_inner, G) = apply.get_binary(app_outer)
1078 :     arity_inner= app_inner.opr.arity
1079 :     #print "outer(2) app_inner",app_inner.opr.name
1080 :     if (arity_inner==1):
1081 :     oexp_inner = unary(app_inner)
1082 :     rhs = G
1083 :     (oty_outer, oexp_tmp) = applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
1084 :     return (oty_outer, oexp_tmp)
1085 :     elif(arity_inner==2):
1086 :     oexp_inner = binary(app_inner)
1087 :     #print "applying binary frst time", oexp_inner
1088 :     rhs = G
1089 :     (oty_outer, oexp_tmp) = applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
1090 :     #print "applying binary second time", oexp_tmp
1091 :     return (oty_outer, oexp_tmp)
1092 :     else:
1093 :     raise Exception ("arity is not supported: "+str(arity_outer))
1094 : cchiw 3939 else:
1095 :     raise Exception ("arity is not supported: "+str(arity_outer))
1096 : cchiw 3915
1097 :     # *************************** evaluate at positions ***************************
1098 :     #evaluate scalar field exp
1099 : cchiw 4308 def eval_d0(pos0, exp):
1100 :     return exp
1101 :    
1102 :    
1103 : cchiw 4158 def eval_d1(pos0, exp):
1104 : cchiw 4230 #print "eval vec d1"
1105 : cchiw 4158 #print "exp:",exp
1106 : cchiw 4230 #print "pos0",pos0
1107 : cchiw 4158 #evaluate field defined by coefficients at position
1108 :     exp = exp.subs(x,pos0)
1109 : cchiw 4230 #print "exp",exp
1110 : cchiw 4158 return exp
1111 :    
1112 : cchiw 3915 def eval_d2(pos0, pos1, exp):
1113 : cchiw 3998 #print "exp:",exp
1114 : cchiw 3915 #evaluate field defined by coefficients at position
1115 : cchiw 4158 exp = exp.subs(x,pos0)
1116 :     exp = exp.subs(y,pos1)
1117 :     return exp
1118 : cchiw 3915
1119 : cchiw 3939 def eval_d3(pos0, pos1, pos2, exp):
1120 :     #evaluate field defined by coefficients at position
1121 : cchiw 4158 exp = exp.subs(x,pos0)
1122 :     exp = exp.subs(y,pos1)
1123 :     exp = exp.subs(z,pos2)
1124 :     return exp
1125 : cchiw 3939
1126 : cchiw 4158 #evaluate vector field [exp]
1127 :     def eval_vec_d1(pos0, vec):
1128 : cchiw 4230 #print "eval vec d1"
1129 : cchiw 4158 rtn = []
1130 :     for v in vec:
1131 :     rtn.append(eval_d1(pos0, v))
1132 :     return rtn
1133 :    
1134 : cchiw 3915 #evaluate vector field [exp,exp]
1135 :     def eval_vec_d2(pos0, pos1, vec):
1136 : cchiw 3939 #print "eval_vec_d2 vec:",vec
1137 : cchiw 3915 rtn = []
1138 :     for v in vec:
1139 :     rtn.append(eval_d2(pos0, pos1, v))
1140 :     return rtn
1141 :    
1142 : cchiw 4158 def eval_ten3_d1(pos0,ten3):
1143 :     rtn = []
1144 :     for i in ten3:
1145 :     for j in i:
1146 :     for v in j:
1147 :     rtn.append(eval_d1(pos0, v))
1148 :     return rtn
1149 :    
1150 :    
1151 :    
1152 : cchiw 3939 #evaluate vector field [exp,exp]
1153 : cchiw 4158 def eval_mat_d1(pos0, mat):
1154 :     #print "eval_vec_d2 vec:",vec
1155 :     rtn = []
1156 :     for i in mat:
1157 :     for v in i:
1158 :     rtn.append(eval_d1(pos0, v))
1159 :     return rtn
1160 :    
1161 :     #evaluate vector field [exp,exp]
1162 : cchiw 3946 def eval_mat_d2(pos0, pos1, mat):
1163 :     #print "eval_vec_d2 vec:",vec
1164 :     rtn = []
1165 : cchiw 4210 #print "eval_mat_d2 mat",mat
1166 : cchiw 3946 for i in mat:
1167 :     for v in i:
1168 :     rtn.append(eval_d2(pos0, pos1, v))
1169 :     return rtn
1170 :    
1171 :     def eval_ten3_d2(pos0, pos1, ten3):
1172 :     #print "eval_vec_d2 vec:",vec
1173 :     rtn = []
1174 :     for i in ten3:
1175 :     for j in i:
1176 :     for v in j:
1177 :     rtn.append(eval_d2(pos0, pos1, v))
1178 :     return rtn
1179 :    
1180 :    
1181 :    
1182 :     #evaluate vector field [exp,exp]
1183 : cchiw 3939 def eval_vec_d3(pos0, pos1, pos2, vec):
1184 :     rtn = []
1185 :     for v in vec:
1186 :     rtn.append(eval_d3(pos0, pos1, pos2, v))
1187 :     return rtn
1188 :    
1189 :    
1190 : cchiw 3946 #evaluate vector field [exp,exp]
1191 :     def eval_mat_d3(pos0, pos1, pos2, mat):
1192 :     #print "eval_vec_d2 vec:",vec
1193 :     rtn = []
1194 :     for i in mat:
1195 :     for v in i:
1196 :     rtn.append(eval_d3(pos0, pos1, pos2, v))
1197 :     return rtn
1198 :    
1199 :     def eval_ten3_d3(pos0, pos1, pos2,ten3):
1200 :     rtn = []
1201 :     for i in ten3:
1202 : cchiw 4411
1203 : cchiw 3946 for j in i:
1204 : cchiw 4411
1205 : cchiw 3946 for v in j:
1206 : cchiw 4411
1207 : cchiw 3946 rtn.append(eval_d3(pos0, pos1, pos2, v))
1208 :     return rtn
1209 :    
1210 :    
1211 : cchiw 4158
1212 :     def iter_d1(k, pos, exp):
1213 :     corr = []
1214 :     for x in pos:
1215 :     val = k(x, exp)
1216 :     corr.append(val)
1217 :     return corr
1218 :    
1219 : cchiw 3939 def iter_d2(k, pos, exp):
1220 : cchiw 3874 corr = []
1221 : cchiw 3939 #print "iter expr:", exp
1222 :     #print "pos", pos
1223 : cchiw 3874 for p in pos:
1224 : cchiw 3939 #print "p", p
1225 : cchiw 3874 x=p[0]
1226 :     y=p[1]
1227 : cchiw 3915 val = k(x,y,exp)
1228 : cchiw 3874 corr.append(val)
1229 :     return corr
1230 : cchiw 3865
1231 : cchiw 3939 def iter_d3(k, pos, exp):
1232 :     corr = []
1233 :     #print "iter exp:", exp
1234 :     for p in pos:
1235 :     x=p[0]
1236 :     y=p[1]
1237 :     z=p[2]
1238 :     val = k(x,y,z, exp)
1239 :     #print "pos: ",x,y,z, " val:", val
1240 :     corr.append(val)
1241 :     return corr
1242 :    
1243 : cchiw 3915 def probeField(otyp1,pos, ortn):
1244 : cchiw 3939 dim = fty.get_dim(otyp1)
1245 : cchiw 4411 #print "output type"+otyp1.name
1246 :     #print "inside probe field ortn", ortn
1247 : cchiw 4456 write_outer("dim"+str(dim))
1248 : cchiw 4308 if (dim==nonefield_dim):
1249 :     #still need to flatten
1250 :     rtn = []
1251 :     if (fty.is_Matrix(otyp1)):
1252 :     for i in ortn:
1253 :     for j in i :
1254 :     rtn.append(j)
1255 :     return [rtn]
1256 :     elif (fty.is_Ten3(otyp1)):
1257 :     for i in ortn:
1258 :     for j in i :
1259 :     for k in j:
1260 :     rtn.append(k)
1261 :     return [rtn]
1262 :    
1263 :     else:
1264 :     return [ortn]
1265 :    
1266 :     elif (dim==1):
1267 : cchiw 3939 def get_k():
1268 :     if (fty.is_ScalarField(otyp1)): # output is a scalar field
1269 : cchiw 4230 #print "s_d1"
1270 : cchiw 4158 return eval_d1
1271 :     elif (fty.is_VectorField(otyp1)):
1272 : cchiw 4230 #print "v_d1"
1273 : cchiw 4158 return eval_vec_d1
1274 :     elif (fty.is_Matrix(otyp1)):
1275 :     return eval_mat_d1
1276 :     elif(fty.is_Ten3(otyp1)):
1277 :     return eval_ten3_d1
1278 :     else:
1279 :     raise "error"+otyp1.name
1280 :     return iter_d1(get_k(), pos, ortn)
1281 :     elif (dim==2):
1282 :     def get_k():
1283 :     if (fty.is_ScalarField(otyp1)): # output is a scalar field
1284 : cchiw 3939 return eval_d2
1285 : cchiw 3946 elif (fty.is_VectorField(otyp1)):
1286 :     return eval_vec_d2
1287 :     elif (fty.is_Matrix(otyp1)):
1288 :     return eval_mat_d2
1289 :     elif(fty.is_Ten3(otyp1)):
1290 :     return eval_ten3_d2
1291 : cchiw 3939 else:
1292 : cchiw 3946 raise "error"+otyp1.name
1293 : cchiw 3939 return iter_d2(get_k(), pos, ortn)
1294 :     elif (dim==3):
1295 :     def get_k():
1296 :     if (fty.is_ScalarField(otyp1)): # output is a scalar field
1297 :     return eval_d3
1298 : cchiw 3946 elif (fty.is_VectorField(otyp1)):
1299 :     return eval_vec_d3
1300 :     elif (fty.is_Matrix(otyp1)):
1301 :     return eval_mat_d3
1302 :     elif(fty.is_Ten3(otyp1)):
1303 :     return eval_ten3_d3
1304 : cchiw 3939 else:
1305 : cchiw 3946 raise "error"+otyp1.name
1306 : cchiw 3939 return iter_d3(get_k(), pos, ortn)
1307 : cchiw 3865 else:
1308 : cchiw 3939 raise "unsupported dimension"
1309 : cchiw 3865
1310 : cchiw 3915 # *************************** main ***************************
1311 :    
1312 :     # operators with scalar field and vector field
1313 :     def eval(app, pos):
1314 : cchiw 3946 #print "evalname",app.name
1315 : cchiw 3939 #apply.toStr(app)
1316 : cchiw 3915 (otyp1, ortn) = sort(app) #apply operations to expressions
1317 : cchiw 4210 #print "ortn",ortn
1318 : cchiw 4188 rtn = probeField(otyp1, pos, ortn) #evaluate expression at positions
1319 : cchiw 4210 #print "rtn", rtn
1320 : cchiw 4188 return rtn

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