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
 [diderot] / branches / ein16 / synth / d2 / test_eval.py

# Annotation of /branches/ein16/synth/d2/test_eval.py

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