30 |
|
|
31 |
def generate(DIM): |
def generate(DIM): |
32 |
# variables |
# variables |
33 |
h=[Symbol('h%s'%i) for i in xrange(DIM) ] |
h=[Symbol('h%s'%i) for i in range(DIM) ] |
34 |
x=[Symbol('x%s'%i) for i in xrange(DIM) ] |
x=[Symbol('x%s'%i) for i in range(DIM) ] |
35 |
# shape functions: labeling differently from finley!!!! |
# shape functions: labeling differently from finley!!!! |
36 |
S=[ 1 ] |
S=[ 1 ] |
37 |
GridOffset=[[]] |
GridOffset=[[]] |
41 |
W = [ 1 ] |
W = [ 1 ] |
42 |
Q_r = [ [] ] |
Q_r = [ [] ] |
43 |
W_r = [ 1 ] |
W_r = [ 1 ] |
44 |
for i in xrange(DIM): |
for i in range(DIM): |
45 |
idx_full.append(i) |
idx_full.append(i) |
46 |
S = [ s * (1-x[i]/h[i]) for s in S] + [ s * x[i]/h[i] for s in S] |
S = [ s * (1-x[i]/h[i]) for s in S] + [ s * x[i]/h[i] for s in S] |
47 |
GridOffset = [ s + [0] for s in GridOffset] + [ s + [1] for s in GridOffset] |
GridOffset = [ s + [0] for s in GridOffset] + [ s + [1] for s in GridOffset] |
49 |
|
|
50 |
# generate quadrature points in element |
# generate quadrature points in element |
51 |
Qnew, Wnew =[], [] |
Qnew, Wnew =[], [] |
52 |
for j in xrange(len(qD1)): |
for j in range(len(qD1)): |
53 |
Qnew += [ q + [qD1[j]*h[i],] for q in Q ] |
Qnew += [ q + [qD1[j]*h[i],] for q in Q ] |
54 |
Wnew += [ w * wD1[j]*h[i] for w in W ] |
Wnew += [ w * wD1[j]*h[i] for w in W ] |
55 |
Q, W=Qnew, Wnew |
Q, W=Qnew, Wnew |
56 |
|
|
57 |
Qnew, Wnew =[], [] |
Qnew, Wnew =[], [] |
58 |
for j in xrange(len(qD1_r)): |
for j in range(len(qD1_r)): |
59 |
Qnew += [ q + [qD1_r[j]*h[i],] for q in Q_r ] |
Qnew += [ q + [qD1_r[j]*h[i],] for q in Q_r ] |
60 |
Wnew += [ w * wD1_r[j]*h[i] for w in W_r ] |
Wnew += [ w * wD1_r[j]*h[i] for w in W_r ] |
61 |
Q_r, W_r=Qnew, Wnew |
Q_r, W_r=Qnew, Wnew |
66 |
Q_r_faces=[] |
Q_r_faces=[] |
67 |
W_r_faces=[] |
W_r_faces=[] |
68 |
idx_faces=[] |
idx_faces=[] |
69 |
for face in xrange(DIM): |
for face in range(DIM): |
70 |
Q_left, W_left = [ [] ], [ 1 ] |
Q_left, W_left = [ [] ], [ 1 ] |
71 |
Q_left_r, W_left_r = [ [] ], [ 1 ] |
Q_left_r, W_left_r = [ [] ], [ 1 ] |
72 |
Q_right, W_right = [ [] ], [ 1 ] |
Q_right, W_right = [ [] ], [ 1 ] |
73 |
Q_right_r, W_right_r = [ [] ], [ 1 ] |
Q_right_r, W_right_r = [ [] ], [ 1 ] |
74 |
idx_left,idx_right=[],[] |
idx_left,idx_right=[],[] |
75 |
for i in xrange(DIM): |
for i in range(DIM): |
76 |
# generate quadrature points in element |
# generate quadrature points in element |
77 |
Q_left_new, W_left_new =[], [] |
Q_left_new, W_left_new =[], [] |
78 |
Q_right_new, W_right_new =[], [] |
Q_right_new, W_right_new =[], [] |
87 |
else: |
else: |
88 |
idx_left.append(i) |
idx_left.append(i) |
89 |
idx_right.append(i) |
idx_right.append(i) |
90 |
for j in xrange(len(qD1)): |
for j in range(len(qD1)): |
91 |
Q_left_new += [ q + [qD1[j]*h[i],] for q in Q_left ] |
Q_left_new += [ q + [qD1[j]*h[i],] for q in Q_left ] |
92 |
W_left_new += [ w * wD1[j]*h[i] for w in W_left ] |
W_left_new += [ w * wD1[j]*h[i] for w in W_left ] |
93 |
Q_right_new += [ q + [qD1[j]*h[i],] for q in Q_right ] |
Q_right_new += [ q + [qD1[j]*h[i],] for q in Q_right ] |
102 |
Q_left_r, W_left_r = [ [] ], [ 1 ] |
Q_left_r, W_left_r = [ [] ], [ 1 ] |
103 |
Q_right, W_right = [ [] ], [ 1 ] |
Q_right, W_right = [ [] ], [ 1 ] |
104 |
Q_right_r, W_right_r = [ [] ], [ 1 ] |
Q_right_r, W_right_r = [ [] ], [ 1 ] |
105 |
for i in xrange(DIM): |
for i in range(DIM): |
106 |
# generate quadrature points in element |
# generate quadrature points in element |
107 |
Q_left_new, W_left_new =[], [] |
Q_left_new, W_left_new =[], [] |
108 |
Q_right_new, W_right_new =[], [] |
Q_right_new, W_right_new =[], [] |
112 |
Q_right_new += [ q + [ h[i], ] for q in Q_right ] |
Q_right_new += [ q + [ h[i], ] for q in Q_right ] |
113 |
W_right_new += [ w * 1. for w in W_right ] |
W_right_new += [ w * 1. for w in W_right ] |
114 |
else: |
else: |
115 |
for j in xrange(len(qD1_r)): |
for j in range(len(qD1_r)): |
116 |
Q_left_new += [ q + [qD1_r[j]*h[i],] for q in Q_left ] |
Q_left_new += [ q + [qD1_r[j]*h[i],] for q in Q_left ] |
117 |
W_left_new += [ w * wD1_r[j]*h[i] for w in W_left ] |
W_left_new += [ w * wD1_r[j]*h[i] for w in W_left ] |
118 |
Q_right_new += [ q + [qD1_r[j]*h[i],] for q in Q_right ] |
Q_right_new += [ q + [qD1_r[j]*h[i],] for q in Q_right ] |
124 |
|
|
125 |
# the interpolation function |
# the interpolation function |
126 |
f_interpolation = 0 |
f_interpolation = 0 |
127 |
for i in xrange(len(GridOffsetText)): |
for i in range(len(GridOffsetText)): |
128 |
f_interpolation = f_interpolation + S[i] * Symbol("f_%s"%GridOffsetText[i]) |
f_interpolation = f_interpolation + S[i] * Symbol("f_%s"%GridOffsetText[i]) |
129 |
|
|
130 |
# interpolation to Quad points |
# interpolation to Quad points |
133 |
CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n" |
CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n" |
134 |
CODE+=createCode(f_interpolation, x, Q_r, loopindex=idx_full) |
CODE+=createCode(f_interpolation, x, Q_r, loopindex=idx_full) |
135 |
CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n" |
CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n" |
136 |
for i in xrange(len(Q_r_faces)): |
for i in range(len(Q_r_faces)): |
137 |
CODE+="if (face_offset(%s)>-1) {\n"%i |
CODE+="if (face_offset(%s)>-1) {\n"%i |
138 |
CODE+=createCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) |
CODE+=createCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) |
139 |
CODE+="\n} /* end of face %s */\n"%i |
CODE+="\n} /* end of face %s */\n"%i |
140 |
CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n" |
CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n" |
141 |
for i in xrange(len(Q_r_faces)): |
for i in range(len(Q_r_faces)): |
142 |
CODE+="if (face_offset(%s)>-1) {\n"%i |
CODE+="if (face_offset(%s)>-1) {\n"%i |
143 |
CODE+=createCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) |
CODE+=createCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) |
144 |
CODE+="\n} /* end of face %s */\n"%i |
CODE+="\n} /* end of face %s */\n"%i |
151 |
CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n" |
CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n" |
152 |
CODE+=createGradientCode(f_interpolation, x, Q_r, loopindex=idx_full) |
CODE+=createGradientCode(f_interpolation, x, Q_r, loopindex=idx_full) |
153 |
CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n" |
CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n" |
154 |
for i in xrange(len(Q_r_faces)): |
for i in range(len(Q_r_faces)): |
155 |
CODE+="if (face_offset(%s)>-1) {\n"%i |
CODE+="if (face_offset(%s)>-1) {\n"%i |
156 |
CODE+=createGradientCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) |
CODE+=createGradientCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) |
157 |
CODE+="\n} /* end of face %s */\n"%i |
CODE+="\n} /* end of face %s */\n"%i |
158 |
CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n" |
CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n" |
159 |
for i in xrange(len(Q_r_faces)): |
for i in range(len(Q_r_faces)): |
160 |
CODE+="if (face_offset(%s)>-1) {\n"%i |
CODE+="if (face_offset(%s)>-1) {\n"%i |
161 |
CODE+=createGradientCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) |
CODE+=createGradientCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) |
162 |
CODE+="\n} /* end of face %s */\n"%i |
CODE+="\n} /* end of face %s */\n"%i |
196 |
BOTTOM_LINE=BOTTOM_LINE_FMT%k |
BOTTOM_LINE=BOTTOM_LINE_FMT%k |
197 |
new_code="" |
new_code="" |
198 |
ignore=False |
ignore=False |
199 |
for l in old_code.split("\n"): |
for l in old_code.split("\n"): |
200 |
s=l.strip() |
s=l.strip() |
201 |
if len(s)>0: |
if len(s)>0: |
202 |
if s.startswith(TOP_LINE): |
if s.startswith(TOP_LINE): |
214 |
for l in old_code.split("\n"): |
for l in old_code.split("\n"): |
215 |
s=l.strip() |
s=l.strip() |
216 |
if s.find("/*")>-1 and s.find("*/")>-1 : |
if s.find("/*")>-1 and s.find("*/")>-1 : |
217 |
s2=(s[:s.find("/*")] + s[s.find("*/")+2:]).strip() |
s2=(s[:s.find("/*")] + s[s.find("*/")+2:]).strip() |
218 |
else: |
else: |
219 |
s2=s |
s2=s |
220 |
if len(s)>0: |
if len(s)>0: |
221 |
if s2.startswith("}"): |
if s2.startswith("}"): |
222 |
N=max(N-1, 0) |
N=max(N-1, 0) |
232 |
F2=[] |
F2=[] |
233 |
for f in F: |
for f in F: |
234 |
f0=[] |
f0=[] |
235 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
236 |
tmp=f |
tmp=f |
237 |
for i in range(len(x)): tmp=tmp.subs(x[i], Q[q][i]) |
for i in range(len(x)): tmp=tmp.subs(x[i], Q[q][i]) |
238 |
f0.append(tmp) |
f0.append(tmp) |
281 |
k="" |
k="" |
282 |
N="" |
N="" |
283 |
M1="" |
M1="" |
284 |
for i in xrange(len(Q[0])): |
for i in range(len(Q[0])): |
285 |
if len(k)==0: |
if len(k)==0: |
286 |
k+="k%s"%i |
k+="k%s"%i |
287 |
N+="N%s"%i |
N+="N%s"%i |
292 |
if i < DIM-1: M1+=",M%s"%i |
if i < DIM-1: M1+=",M%s"%i |
293 |
|
|
294 |
k3="" |
k3="" |
295 |
for i in xrange(DIM-1,-1,-1): |
for i in range(DIM-1,-1,-1): |
296 |
if loopindex[i] > -1: k3+=",k%i"%loopindex[i] |
if loopindex[i] > -1: k3+=",k%i"%loopindex[i] |
297 |
TXT="" |
TXT="" |
298 |
for a,v in consts.items(): |
for a,v in consts.items(): |
299 |
TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS))) |
TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS))) |
300 |
TXT+="#pragma omp parallel for private(i%s)"%k3 |
TXT+="#pragma omp parallel for private(i%s)"%k3 |
301 |
|
|
302 |
for i in xrange(DIM-1,-1,-1): |
for i in range(DIM-1,-1,-1): |
303 |
if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i]) |
if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i]) |
304 |
TXT+="\nfor (i =0; i < NCOMP; ++i) {\n" |
TXT+="\nfor (i =0; i < NCOMP; ++i) {\n" |
305 |
for s in args: |
for s in args: |
306 |
k1="" |
k1="" |
307 |
for i in xrange(DIM): |
for i in range(DIM): |
308 |
if s.name[2+i]=="0": |
if s.name[2+i]=="0": |
309 |
if loopindex[i]>-1: |
if loopindex[i]>-1: |
310 |
k1+=",k%s"%i |
k1+=",k%s"%i |
321 |
k1+=",M%s-1"%i |
k1+=",M%s-1"%i |
322 |
TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1) |
TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1) |
323 |
# interpolation to quadrature points |
# interpolation to quadrature points |
324 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
325 |
IDX2="INDEX%s(%s,%s)"%(DIM,k,N) |
IDX2="INDEX%s(%s,%s)"%(DIM,k,N) |
326 |
if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2 |
if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2 |
327 |
for i in range(DIM): |
for i in range(DIM): |
328 |
TXT+="out[INDEX4(i,%s,%s,%s,NCOMP,%s,%s)] = %s;\n"%(i,q,IDX2,DIM,len(Q),ccode(dF[i][q])) |
TXT+="out[INDEX4(i,%s,%s,%s,NCOMP,%s,%s)] = %s;\n"%(i,q,IDX2,DIM,len(Q),ccode(dF[i][q])) |
329 |
TXT+="} /* close component loop i */\n" |
TXT+="} /* close component loop i */\n" |
330 |
for i in xrange(DIM): |
for i in range(DIM): |
331 |
if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i] |
if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i] |
332 |
return TXT |
return TXT |
333 |
|
|
338 |
k="" |
k="" |
339 |
N="" |
N="" |
340 |
M1="" |
M1="" |
341 |
for i in xrange(len(Q[0])): |
for i in range(len(Q[0])): |
342 |
if len(k)==0: |
if len(k)==0: |
343 |
k+="k%s"%i |
k+="k%s"%i |
344 |
N+="N%s"%i |
N+="N%s"%i |
349 |
if i < DIM-1: M1+=",M%s"%i |
if i < DIM-1: M1+=",M%s"%i |
350 |
|
|
351 |
k3="" |
k3="" |
352 |
for i in xrange(DIM-1,-1,-1): |
for i in range(DIM-1,-1,-1): |
353 |
if loopindex[i] > -1: k3+=",k%i"%loopindex[i] |
if loopindex[i] > -1: k3+=",k%i"%loopindex[i] |
354 |
TXT="" |
TXT="" |
355 |
for a,v in consts.items(): |
for a,v in consts.items(): |
356 |
TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS))) |
TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS))) |
357 |
TXT+="#pragma omp parallel for private(i%s)"%k3 |
TXT+="#pragma omp parallel for private(i%s)"%k3 |
358 |
for i in xrange(DIM-1,-1,-1): |
for i in range(DIM-1,-1,-1): |
359 |
if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i]) |
if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i]) |
360 |
TXT+="\nfor (i =0; i < NCOMP; ++i) {\n" |
TXT+="\nfor (i =0; i < NCOMP; ++i) {\n" |
361 |
for s in args: |
for s in args: |
362 |
k1="" |
k1="" |
363 |
for i in xrange(DIM): |
for i in range(DIM): |
364 |
if s.name[2+i]=="0": |
if s.name[2+i]=="0": |
365 |
if loopindex[i]>-1: |
if loopindex[i]>-1: |
366 |
k1+=",k%s"%i |
k1+=",k%s"%i |
377 |
k1+=",M%s-1"%i |
k1+=",M%s-1"%i |
378 |
TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1) |
TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1) |
379 |
# interpolation to quadrature points |
# interpolation to quadrature points |
380 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
381 |
IDX2="INDEX%s(%s,%s)"%(DIM,k,N) |
IDX2="INDEX%s(%s,%s)"%(DIM,k,N) |
382 |
if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2 |
if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2 |
383 |
TXT+="out[INDEX3(i,%s,%s,NCOMP,%s)] = %s;\n"%(q,IDX2,len(Q),ccode(F[q])) |
TXT+="out[INDEX3(i,%s,%s,NCOMP,%s)] = %s;\n"%(q,IDX2,len(Q),ccode(F[q])) |
384 |
TXT+="} /* close component loop i */\n" |
TXT+="} /* close component loop i */\n" |
385 |
for i in xrange(DIM): |
for i in range(DIM): |
386 |
if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i] |
if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i] |
387 |
return TXT |
return TXT |
388 |
|
|
392 |
# return out.evalf() |
# return out.evalf() |
393 |
return out |
return out |
394 |
|
|
395 |
def generatePDECode(DATA_A, EM,GLOBAL_TMP, system=False): |
def generatePDECode(DATA_A, EM,GLOBAL_TMP, system=False): |
396 |
LOCAL_TMP={} |
LOCAL_TMP={} |
397 |
LOCAL2_TMP={} |
LOCAL2_TMP={} |
398 |
OUT="" |
OUT="" |
399 |
for p in EM.keys(): |
for p in EM.keys(): |
400 |
# first we collect the constant coefficients in the expression (outside the parallel region): |
# first we collect the constant coefficients in the expression (outside the parallel region): |
401 |
tmp={} |
tmp={} |
402 |
for a in DATA_A: |
for a in DATA_A: |
403 |
c= collect(EM[p], a, evaluate=False) |
c= collect(EM[p], a, evaluate=False) |
404 |
if c.has_key(a): |
if c.has_key(a): |
405 |
if not GLOBAL_TMP.has_key(c[a]): |
if not GLOBAL_TMP.has_key(c[a]): |
406 |
GLOBAL_TMP[c[a]]=Symbol("w%s"%(len(GLOBAL_TMP),)) |
GLOBAL_TMP[c[a]]=Symbol("w%s"%(len(GLOBAL_TMP),)) |
407 |
if tmp.has_key(GLOBAL_TMP[c[a]]): |
if tmp.has_key(GLOBAL_TMP[c[a]]): |
408 |
tmp[GLOBAL_TMP[c[a]]]+=a |
tmp[GLOBAL_TMP[c[a]]]+=a |
409 |
else: |
else: |
410 |
tmp[GLOBAL_TMP[c[a]]]=a |
tmp[GLOBAL_TMP[c[a]]]=a |
411 |
# now we find temporary values which are combinations of the input data: |
# now we find temporary values which are combinations of the input data: |
412 |
em_new=[] |
em_new=[] |
413 |
for t in tmp.keys(): |
for t in tmp.keys(): |
414 |
if tmp[t] in DATA_A: |
if tmp[t] in DATA_A: |
415 |
tt = t * tmp[t] |
tt = t * tmp[t] |
416 |
else: |
else: |
417 |
if not LOCAL_TMP.has_key(tmp[t]): |
if not LOCAL_TMP.has_key(tmp[t]): |
418 |
LOCAL_TMP[tmp[t]]=Symbol("tmp%s_0"%len(LOCAL_TMP)) |
LOCAL_TMP[tmp[t]]=Symbol("tmp%s_0"%len(LOCAL_TMP)) |
419 |
tt= t * LOCAL_TMP[tmp[t]] |
tt= t * LOCAL_TMP[tmp[t]] |
420 |
if not LOCAL2_TMP.has_key(tt): |
if not LOCAL2_TMP.has_key(tt): |
421 |
LOCAL2_TMP[tt]=Symbol("tmp%s_1"%len(LOCAL2_TMP)) |
LOCAL2_TMP[tt]=Symbol("tmp%s_1"%len(LOCAL2_TMP)) |
422 |
em_new.append(LOCAL2_TMP[tt]) |
em_new.append(LOCAL2_TMP[tt]) |
423 |
EM[p]=em_new |
EM[p]=em_new |
424 |
|
|
425 |
for p in EM.keys(): |
for p in EM.keys(): |
426 |
sum_em=0 |
sum_em=0 |
427 |
for t in EM[p]: sum_em=sum_em + t |
for t in EM[p]: sum_em=sum_em + t |
428 |
if isinstance(p, int) : |
if isinstance(p, int) : |
429 |
if system: |
if system: |
430 |
OUT+=" EM_F[INDEX2(k,%s,p.numEqu)]+=%s;\n"%(p,ccode(sum_em)) |
OUT+=" EM_F[INDEX2(k,%s,p.numEqu)]+=%s;\n"%(p,ccode(sum_em)) |
431 |
else: |
else: |
432 |
OUT+=" EM_F[%s]+=%s;\n"%(p,ccode(sum_em)) |
OUT+=" EM_F[%s]+=%s;\n"%(p,ccode(sum_em)) |
433 |
else: |
else: |
434 |
if system: |
if system: |
435 |
OUT+=" EM_S[INDEX4(k,m,%s,%s,p.numEqu,p.numComp,%s)]+=%s;\n"%(p[0],p[1],max([ qq[0] for qq in EM.keys()])+1,ccode(sum_em)) |
OUT+=" EM_S[INDEX4(k,m,%s,%s,p.numEqu,p.numComp,%s)]+=%s;\n"%(p[0],p[1],max([ qq[0] for qq in EM.keys()])+1,ccode(sum_em)) |
436 |
else: |
else: |
437 |
OUT+=" EM_S[INDEX2(%s,%s,%s)]+=%s;\n"%(p[0],p[1],max([ qq[0] for qq in EM.keys()])+1,ccode(sum_em)) |
OUT+=" EM_S[INDEX2(%s,%s,%s)]+=%s;\n"%(p[0],p[1],max([ qq[0] for qq in EM.keys()])+1,ccode(sum_em)) |
438 |
|
|
439 |
OUT2="" |
OUT2="" |
440 |
for p in LOCAL_TMP: |
for p in LOCAL_TMP: |
441 |
OUT2+=" const register double %s = %s;\n"%(LOCAL_TMP[p],ccode(p)) |
OUT2+=" const register double %s = %s;\n"%(LOCAL_TMP[p],ccode(p)) |
442 |
for p in LOCAL2_TMP: |
for p in LOCAL2_TMP: |
443 |
OUT2+=" const register double %s = %s;\n"%(LOCAL2_TMP[p],ccode(p)) |
OUT2+=" const register double %s = %s;\n"%(LOCAL2_TMP[p],ccode(p)) |
444 |
return OUT2+OUT |
return OUT2+OUT |
445 |
|
|
446 |
def makePDE(S, x, Q, W, DIM=2, system=False): |
def makePDE(S, x, Q, W, DIM=2, system=False): |
463 |
CODE2="" |
CODE2="" |
464 |
EM = {} |
EM = {} |
465 |
for i in range(len(S)): |
for i in range(len(S)): |
466 |
for j in range(len(S)): |
for j in range(len(S)): |
467 |
EM[(i,j)] = 0 |
EM[(i,j)] = 0 |
468 |
|
|
469 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
470 |
for di in range(DIM): |
for di in range(DIM): |
471 |
for dj in range(DIM): |
for dj in range(DIM): |
472 |
A_name="A_%d%d_%d"%(di,dj,q) |
A_name="A_%d%d_%d"%(di,dj,q) |
473 |
A=Symbol(A_name) |
A=Symbol(A_name) |
474 |
DATA_A.append(A) |
DATA_A.append(A) |
475 |
if system: |
if system: |
476 |
CODE2+="const register double %s = A_p[INDEX5(k,%s,m,%s,%s, p.numEqu,%s,p.numComp,%s)];\n"%(A_name, di, dj, q, DIM, DIM) |
CODE2+="const register double %s = A_p[INDEX5(k,%s,m,%s,%s, p.numEqu,%s,p.numComp,%s)];\n"%(A_name, di, dj, q, DIM, DIM) |
477 |
else: |
else: |
478 |
CODE2+="const register double %s = A_p[INDEX3(%s,%s,%s,%s,%s)];\n"%(A_name, di, dj, q, DIM, DIM) |
CODE2+="const register double %s = A_p[INDEX3(%s,%s,%s,%s,%s)];\n"%(A_name, di, dj, q, DIM, DIM) |
479 |
for i in range(len(S)): |
for i in range(len(S)): |
480 |
for j in range(len(S)): |
for j in range(len(S)): |
481 |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[i],x[di]) * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[i],x[di]) * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
482 |
CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system) |
483 |
|
|
484 |
if system: CODE+="}\n }\n" |
if system: CODE+="}\n }\n" |
485 |
CODE+="} else { /* constant data */\n" |
CODE+="} else { /* constant data */\n" |
486 |
if system: |
if system: |
487 |
CODE+= """for (k=0;k<p.numEqu;k++) { |
CODE+= """for (k=0;k<p.numEqu;k++) { |
488 |
for (m=0;m<p.numComp;m++) { |
for (m=0;m<p.numComp;m++) { |
492 |
EM = {} |
EM = {} |
493 |
for i in range(len(S)): |
for i in range(len(S)): |
494 |
for j in range(len(S)): |
for j in range(len(S)): |
495 |
EM[(i,j)] = 0 |
EM[(i,j)] = 0 |
496 |
for di in range(DIM): |
for di in range(DIM): |
497 |
for dj in range(DIM): |
for dj in range(DIM): |
498 |
A_name="A_%d%d"%(di,dj) |
A_name="A_%d%d"%(di,dj) |
499 |
A=Symbol(A_name) |
A=Symbol(A_name) |
500 |
DATA_A.append(A) |
DATA_A.append(A) |
501 |
if system: |
if system: |
502 |
CODE2+="const register double %s = A_p[INDEX4(k,%s,m,%s p.numEqu,%s, p.numComp)];\n"%(A_name, di, dj, DIM) |
CODE2+="const register double %s = A_p[INDEX4(k,%s,m,%s p.numEqu,%s, p.numComp)];\n"%(A_name, di, dj, DIM) |
503 |
else: |
else: |
504 |
CODE2+="const register double %s = A_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, dj, DIM) |
CODE2+="const register double %s = A_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, dj, DIM) |
505 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
506 |
for i in range(len(S)): |
for i in range(len(S)): |
507 |
for j in range(len(S)): |
for j in range(len(S)): |
508 |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[i],x[di]) * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[i],x[di]) * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
509 |
CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system) |
510 |
if system: CODE+="}\n }\n" |
if system: CODE+="}\n }\n" |
511 |
if len(Q) > 1: CODE+="}\n" |
if len(Q) > 1: CODE+="}\n" |
527 |
CODE2="" |
CODE2="" |
528 |
EM = {} |
EM = {} |
529 |
for i in range(len(S)): |
for i in range(len(S)): |
530 |
for j in range(len(S)): |
for j in range(len(S)): |
531 |
EM[(i,j)] = 0 |
EM[(i,j)] = 0 |
532 |
|
|
533 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
534 |
for di in range(DIM): |
for di in range(DIM): |
535 |
A_name="B_%d_%d"%(di,q) |
A_name="B_%d_%d"%(di,q) |
536 |
A=Symbol(A_name) |
A=Symbol(A_name) |
537 |
DATA_B.append(A) |
DATA_B.append(A) |
538 |
if system: |
if system: |
539 |
CODE2+="const register double %s = B_p[INDEX4(k,%s,m,%s, p.numEqu,%s,p.numComp)];\n"%(A_name, di, q, DIM) |
CODE2+="const register double %s = B_p[INDEX4(k,%s,m,%s, p.numEqu,%s,p.numComp)];\n"%(A_name, di, q, DIM) |
540 |
else: |
else: |
541 |
CODE2+="const register double %s = B_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, q, DIM) |
CODE2+="const register double %s = B_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, q, DIM) |
542 |
for i in range(len(S)): |
for i in range(len(S)): |
543 |
for j in range(len(S)): |
for j in range(len(S)): |
544 |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[i],x[di]) * S[j]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[i],x[di]) * S[j]).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
545 |
CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system) |
546 |
|
|
547 |
if system: CODE+="}\n }\n" |
if system: CODE+="}\n }\n" |
548 |
CODE+="} else { /* constant data */\n" |
CODE+="} else { /* constant data */\n" |
549 |
if system: |
if system: |
550 |
CODE+= """for (k=0;k<p.numEqu;k++) { |
CODE+= """for (k=0;k<p.numEqu;k++) { |
551 |
for (m=0;m<p.numComp;m++) { |
for (m=0;m<p.numComp;m++) { |
555 |
EM = {} |
EM = {} |
556 |
for i in range(len(S)): |
for i in range(len(S)): |
557 |
for j in range(len(S)): |
for j in range(len(S)): |
558 |
EM[(i,j)] = 0 |
EM[(i,j)] = 0 |
559 |
for di in range(DIM): |
for di in range(DIM): |
560 |
A_name="B_%d"%(di) |
A_name="B_%d"%(di) |
561 |
A=Symbol(A_name) |
A=Symbol(A_name) |
562 |
DATA_B.append(A) |
DATA_B.append(A) |
563 |
if system: |
if system: |
564 |
CODE2+="const register double %s = B_p[INDEX3(k,%s,m, p.numEqu,%s)];\n"%(A_name, di, DIM) |
CODE2+="const register double %s = B_p[INDEX3(k,%s,m, p.numEqu,%s)];\n"%(A_name, di, DIM) |
565 |
else: |
else: |
566 |
CODE2+="const register double %s = B_p[%s];\n"%(A_name, di) |
CODE2+="const register double %s = B_p[%s];\n"%(A_name, di) |
567 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
568 |
for i in range(len(S)): |
for i in range(len(S)): |
569 |
for j in range(len(S)): |
for j in range(len(S)): |
570 |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[i],x[di]) * S[j] ).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[i],x[di]) * S[j] ).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
571 |
CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system) |
572 |
if system: CODE+="}\n }\n" |
if system: CODE+="}\n }\n" |
573 |
if len(Q) > 1: CODE+="}\n" |
if len(Q) > 1: CODE+="}\n" |
589 |
CODE2="" |
CODE2="" |
590 |
EM = {} |
EM = {} |
591 |
for i in range(len(S)): |
for i in range(len(S)): |
592 |
for j in range(len(S)): |
for j in range(len(S)): |
593 |
EM[(i,j)] = 0 |
EM[(i,j)] = 0 |
594 |
|
|
595 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
596 |
for dj in range(DIM): |
for dj in range(DIM): |
597 |
A_name="C_%d_%d"%(dj,q) |
A_name="C_%d_%d"%(dj,q) |
598 |
A=Symbol(A_name) |
A=Symbol(A_name) |
599 |
DATA_C.append(A) |
DATA_C.append(A) |
600 |
if system: |
if system: |
601 |
CODE2+="const register double %s = C_p[INDEX4(k,m,%s, %s, p.numEqu,p.numComp,%s)];\n"%(A_name, dj, q, DIM) |
CODE2+="const register double %s = C_p[INDEX4(k,m,%s, %s, p.numEqu,p.numComp,%s)];\n"%(A_name, dj, q, DIM) |
602 |
else: |
else: |
603 |
CODE2+="const register double %s = C_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj, q, DIM) |
CODE2+="const register double %s = C_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj, q, DIM) |
604 |
for i in range(len(S)): |
for i in range(len(S)): |
605 |
for j in range(len(S)): |
for j in range(len(S)): |
606 |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[j],x[dj]) * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[j],x[dj]) * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
607 |
CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system) |
608 |
|
|
609 |
if system: CODE+="}\n }\n" |
if system: CODE+="}\n }\n" |
610 |
CODE+="} else { /* constant data */\n" |
CODE+="} else { /* constant data */\n" |
611 |
if system: |
if system: |
612 |
CODE+= """for (k=0;k<p.numEqu;k++) { |
CODE+= """for (k=0;k<p.numEqu;k++) { |
613 |
for (m=0;m<p.numComp;m++) { |
for (m=0;m<p.numComp;m++) { |
617 |
EM = {} |
EM = {} |
618 |
for i in range(len(S)): |
for i in range(len(S)): |
619 |
for j in range(len(S)): |
for j in range(len(S)): |
620 |
EM[(i,j)] = 0 |
EM[(i,j)] = 0 |
621 |
for dj in range(DIM): |
for dj in range(DIM): |
622 |
A_name="C_%d"%(dj) |
A_name="C_%d"%(dj) |
623 |
A=Symbol(A_name) |
A=Symbol(A_name) |
624 |
DATA_C.append(A) |
DATA_C.append(A) |
625 |
if system: |
if system: |
626 |
CODE2+="const register double %s = C_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, dj) |
CODE2+="const register double %s = C_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, dj) |
627 |
else: |
else: |
628 |
CODE2+="const register double %s = C_p[%s];\n"%(A_name, dj) |
CODE2+="const register double %s = C_p[%s];\n"%(A_name, dj) |
629 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
630 |
for i in range(len(S)): |
for i in range(len(S)): |
631 |
for j in range(len(S)): |
for j in range(len(S)): |
632 |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[j],x[dj]) * S[i] ).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * diff(S[j],x[dj]) * S[i] ).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
633 |
CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system) |
634 |
if system: CODE+="}\n }\n" |
if system: CODE+="}\n }\n" |
635 |
if len(Q) > 1: CODE+="}\n" |
if len(Q) > 1: CODE+="}\n" |
651 |
CODE2="" |
CODE2="" |
652 |
EM = {} |
EM = {} |
653 |
for i in range(len(S)): |
for i in range(len(S)): |
654 |
for j in range(len(S)): |
for j in range(len(S)): |
655 |
EM[(i,j)] = 0 |
EM[(i,j)] = 0 |
656 |
|
|
657 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
658 |
A_name="D_%d"%(q,) |
A_name="D_%d"%(q,) |
659 |
A=Symbol(A_name) |
A=Symbol(A_name) |
660 |
DATA_D.append(A) |
DATA_D.append(A) |
661 |
if system: |
if system: |
662 |
CODE2+="const register double %s = D_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, q) |
CODE2+="const register double %s = D_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, q) |
663 |
else: |
else: |
664 |
CODE2+="const register double %s = D_p[%s];\n"%(A_name, q) |
CODE2+="const register double %s = D_p[%s];\n"%(A_name, q) |
665 |
for i in range(len(S)): |
for i in range(len(S)): |
666 |
for j in range(len(S)): |
for j in range(len(S)): |
667 |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
668 |
CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system) |
669 |
|
|
670 |
if system: CODE+="}\n }\n" |
if system: CODE+="}\n }\n" |
671 |
CODE+="} else { /* constant data */\n" |
CODE+="} else { /* constant data */\n" |
672 |
if system: |
if system: |
673 |
CODE+= """for (k=0;k<p.numEqu;k++) { |
CODE+= """for (k=0;k<p.numEqu;k++) { |
674 |
for (m=0;m<p.numComp;m++) { |
for (m=0;m<p.numComp;m++) { |
678 |
EM = {} |
EM = {} |
679 |
for i in range(len(S)): |
for i in range(len(S)): |
680 |
for j in range(len(S)): |
for j in range(len(S)): |
681 |
EM[(i,j)] = 0 |
EM[(i,j)] = 0 |
682 |
if 1: |
if 1: |
683 |
A_name="D_0" |
A_name="D_0" |
684 |
A=Symbol(A_name) |
A=Symbol(A_name) |
685 |
DATA_D.append(A) |
DATA_D.append(A) |
686 |
if system: |
if system: |
687 |
CODE2+="const register double %s = D_p[INDEX2(k,m, p.numEqu)];\n"%(A_name,) |
CODE2+="const register double %s = D_p[INDEX2(k,m, p.numEqu)];\n"%(A_name,) |
688 |
else: |
else: |
689 |
CODE2+="const register double %s = D_p[0];\n"%(A_name,) |
CODE2+="const register double %s = D_p[0];\n"%(A_name,) |
690 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
691 |
for i in range(len(S)): |
for i in range(len(S)): |
692 |
for j in range(len(S)): |
for j in range(len(S)): |
693 |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i] ).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i] ).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
694 |
CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system) |
695 |
if system: CODE+="}\n }\n" |
if system: CODE+="}\n }\n" |
696 |
if len(Q) > 1: CODE+="}\n" |
if len(Q) > 1: CODE+="}\n" |
712 |
CODE2="" |
CODE2="" |
713 |
EM = {} |
EM = {} |
714 |
for i in range(len(S)): |
for i in range(len(S)): |
715 |
EM[i] = 0 |
EM[i] = 0 |
716 |
|
|
717 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
718 |
for dj in range(DIM): |
for dj in range(DIM): |
719 |
A_name="X_%d_%d"%(dj,q) |
A_name="X_%d_%d"%(dj,q) |
720 |
A=Symbol(A_name) |
A=Symbol(A_name) |
721 |
DATA_X.append(A) |
DATA_X.append(A) |
722 |
if system: |
if system: |
723 |
CODE2+="const register double %s = X_p[INDEX3(k,%s, %s, p.numEqu,%s)];\n"%(A_name, dj, q, DIM) |
CODE2+="const register double %s = X_p[INDEX3(k,%s, %s, p.numEqu,%s)];\n"%(A_name, dj, q, DIM) |
724 |
else: |
else: |
725 |
CODE2+="const register double %s = X_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj,q,DIM) |
CODE2+="const register double %s = X_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj,q,DIM) |
726 |
for j in range(len(S)): |
for j in range(len(S)): |
727 |
EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
728 |
CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system) |
729 |
|
|
730 |
if system: CODE+="}\n" |
if system: CODE+="}\n" |
731 |
CODE+="} else { /* constant data */\n" |
CODE+="} else { /* constant data */\n" |
732 |
if system: |
if system: |
733 |
CODE+= """for (k=0;k<p.numEqu;k++) { |
CODE+= """for (k=0;k<p.numEqu;k++) { |
734 |
""" |
""" |
736 |
CODE2="" |
CODE2="" |
737 |
EM = {} |
EM = {} |
738 |
for i in range(len(S)): |
for i in range(len(S)): |
739 |
EM[i] = 0 |
EM[i] = 0 |
740 |
for dj in range(DIM): |
for dj in range(DIM): |
741 |
A_name="X_%d"%(dj) |
A_name="X_%d"%(dj) |
742 |
A=Symbol(A_name) |
A=Symbol(A_name) |
743 |
DATA_X.append(A) |
DATA_X.append(A) |
744 |
if system: |
if system: |
745 |
CODE2+="const register double %s = X_p[INDEX2(k,%s, p.numEqu)];\n"%(A_name, dj) |
CODE2+="const register double %s = X_p[INDEX2(k,%s, p.numEqu)];\n"%(A_name, dj) |
746 |
else: |
else: |
747 |
CODE2+="const register double %s = X_p[%s];\n"%(A_name, dj) |
CODE2+="const register double %s = X_p[%s];\n"%(A_name, dj) |
748 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
749 |
for j in range(len(S)): |
for j in range(len(S)): |
750 |
EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
751 |
CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system) |
752 |
if system: CODE+="}\n" |
if system: CODE+="}\n" |
753 |
if len(Q) > 1: CODE+="}\n" |
if len(Q) > 1: CODE+="}\n" |
769 |
CODE2="" |
CODE2="" |
770 |
EM = {} |
EM = {} |
771 |
for i in range(len(S)): |
for i in range(len(S)): |
772 |
EM[i] = 0 |
EM[i] = 0 |
773 |
|
|
774 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
775 |
A_name="Y_%d"%(q,) |
A_name="Y_%d"%(q,) |
776 |
A=Symbol(A_name) |
A=Symbol(A_name) |
777 |
DATA_Y.append(A) |
DATA_Y.append(A) |
778 |
if system: |
if system: |
779 |
CODE2+="const register double %s = Y_p[INDEX3(k,%s, p.numEqu)];\n"%(A_name, q) |
CODE2+="const register double %s = Y_p[INDEX3(k,%s, p.numEqu)];\n"%(A_name, q) |
780 |
else: |
else: |
781 |
CODE2+="const register double %s = Y_p[%s];\n"%(A_name, q) |
CODE2+="const register double %s = Y_p[%s];\n"%(A_name, q) |
782 |
for i in range(len(S)): |
for i in range(len(S)): |
783 |
EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
784 |
CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system) |
785 |
|
|
786 |
if system: CODE+="}\n" |
if system: CODE+="}\n" |
787 |
CODE+="} else { /* constant data */\n" |
CODE+="} else { /* constant data */\n" |
788 |
if system: |
if system: |
789 |
CODE+= """for (k=0;k<p.numEqu;k++) { |
CODE+= """for (k=0;k<p.numEqu;k++) { |
790 |
""" |
""" |
792 |
CODE2="" |
CODE2="" |
793 |
EM = {} |
EM = {} |
794 |
for i in range(len(S)): |
for i in range(len(S)): |
795 |
EM[i] = 0 |
EM[i] = 0 |
796 |
if 1: |
if 1: |
797 |
A_name="Y_0" |
A_name="Y_0" |
798 |
A=Symbol(A_name) |
A=Symbol(A_name) |
799 |
DATA_Y.append(A) |
DATA_Y.append(A) |
800 |
if system: |
if system: |
801 |
CODE2+="const register double %s = Y_p[k];\n"%(A_name,) |
CODE2+="const register double %s = Y_p[k];\n"%(A_name,) |
802 |
else: |
else: |
803 |
CODE2+="const register double %s = Y_p[0];\n"%(A_name,) |
CODE2+="const register double %s = Y_p[0];\n"%(A_name,) |
804 |
for q in xrange(len(Q)): |
for q in range(len(Q)): |
805 |
for i in range(len(S)): |
for i in range(len(S)): |
806 |
EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] ) |
EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] ) |
807 |
CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system) |
CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system) |
808 |
if system: CODE+="}\n" |
if system: CODE+="}\n" |
809 |
if len(Q) > 1: CODE+="}\n" |
if len(Q) > 1: CODE+="}\n" |