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

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