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

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