/[escript]/trunk/ripley/src/generate_assamblage.py
ViewVC logotype

Contents of /trunk/ripley/src/generate_assamblage.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3792 - (show annotations)
Wed Feb 1 06:16:25 2012 UTC (6 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 29666 byte(s)
Merged ripley rectangular domain into trunk.

1 """
2 this script generates the assemblage routine for the ripley rectangular grid solver
3
4
5 """
6 from sympy import *
7
8 DIGITS=20
9
10 #
11 # quadrature form in [0,1]
12 #
13
14 q0=1/sqrt(RealNumber(3))
15 qD1=[ (1-q0)/2, (1+q0)/2 ]
16 wD1=[RealNumber(1)/2, RealNumber(1)/2 ]
17
18
19 q0=1/sqrt(RealNumber(3))
20 qD1=[ (1-q0)/2, (1+q0)/2 ]
21 wD1=[RealNumber(1)/2, RealNumber(1)/2 ]
22
23 qD1_r=[ RealNumber(1)/2 ]
24 wD1_r= [RealNumber(1)]
25
26 print "1D quadrature nodes =",qD1
27 print "1D quadrature weights =",wD1
28 print "1D reduced quadrature nodes =",qD1_r
29 print "1D reduced quadrature weights =",wD1_r
30
31 def generate(DIM):
32 # variables
33 h=[Symbol('h%s'%i) for i in xrange(DIM) ]
34 x=[Symbol('x%s'%i) for i in xrange(DIM) ]
35 # shape functions: labeling differently from finley!!!!
36 S=[ 1 ]
37 GridOffset=[[]]
38 GridOffsetText=[ "" ]
39 idx_full=[]
40 Q = [ [] ]
41 W = [ 1 ]
42 Q_r = [ [] ]
43 W_r = [ 1 ]
44 for i in xrange(DIM):
45 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]
47 GridOffset = [ s + [0] for s in GridOffset] + [ s + [1] for s in GridOffset]
48 GridOffsetText = [ s + "0" for s in GridOffsetText] + [ s + "1" for s in GridOffsetText]
49
50 # generate quadrature points in element
51 Qnew, Wnew =[], []
52 for j in xrange(len(qD1)):
53 Qnew += [ q + [qD1[j]*h[i],] for q in Q ]
54 Wnew += [ w * wD1[j]*h[i] for w in W ]
55 Q, W=Qnew, Wnew
56
57 Qnew, Wnew =[], []
58 for j in xrange(len(qD1_r)):
59 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 ]
61 Q_r, W_r=Qnew, Wnew
62
63 # now the same thing for the faces:
64 Q_faces=[]
65 W_faces=[]
66 Q_r_faces=[]
67 W_r_faces=[]
68 idx_faces=[]
69 for face in xrange(DIM):
70 Q_left, W_left = [ [] ], [ 1 ]
71 Q_left_r, W_left_r = [ [] ], [ 1 ]
72 Q_right, W_right = [ [] ], [ 1 ]
73 Q_right_r, W_right_r = [ [] ], [ 1 ]
74 idx_left,idx_right=[],[]
75 for i in xrange(DIM):
76 # generate quadrature points in element
77 Q_left_new, W_left_new =[], []
78 Q_right_new, W_right_new =[], []
79
80 if face == i :
81 idx_left.append(-1)
82 idx_right.append(-2)
83 Q_left_new += [ q + [0.,] for q in Q_left ]
84 W_left_new += [ w * 1. for w in W_left ]
85 Q_right_new += [ q + [ h[i], ] for q in Q_right ]
86 W_right_new += [ w * 1. for w in W_right ]
87 else:
88 idx_left.append(i)
89 idx_right.append(i)
90 for j in xrange(len(qD1)):
91 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 ]
93 Q_right_new += [ q + [qD1[j]*h[i],] for q in Q_right ]
94 W_right_new += [ w * wD1[j]*h[i] for w in W_right ]
95 Q_left, W_left=Q_left_new, W_left_new
96 Q_right, W_right=Q_right_new, W_right_new
97 Q_faces=Q_faces+[Q_left, Q_right]
98 W_faces=W_faces+[W_left, W_right]
99 idx_faces=idx_faces + [ idx_left, idx_right ]
100
101 Q_left, W_left = [ [] ], [ 1 ]
102 Q_left_r, W_left_r = [ [] ], [ 1 ]
103 Q_right, W_right = [ [] ], [ 1 ]
104 Q_right_r, W_right_r = [ [] ], [ 1 ]
105 for i in xrange(DIM):
106 # generate quadrature points in element
107 Q_left_new, W_left_new =[], []
108 Q_right_new, W_right_new =[], []
109 if face == i :
110 Q_left_new += [ q + [0.,] for q in Q_left ]
111 W_left_new += [ w * 1. for w in W_left ]
112 Q_right_new += [ q + [ h[i], ] for q in Q_right ]
113 W_right_new += [ w * 1. for w in W_right ]
114 else:
115 for j in xrange(len(qD1_r)):
116 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 ]
118 Q_right_new += [ q + [qD1_r[j]*h[i],] for q in Q_right ]
119 W_right_new += [ w * wD1_r[j]*h[i] for w in W_right ]
120 Q_left, W_left=Q_left_new, W_left_new
121 Q_right, W_right=Q_right_new, W_right_new
122 Q_r_faces=Q_r_faces+[Q_left, Q_right]
123 W_r_faces=W_r_faces+[W_left, W_right]
124
125 # the interpolation function
126 f_interpolation = 0
127 for i in xrange(len(GridOffsetText)):
128 f_interpolation = f_interpolation + S[i] * Symbol("f_%s"%GridOffsetText[i])
129
130 # interpolation to Quad points
131 CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"
132 CODE+=createCode(f_interpolation, x, Q, loopindex=idx_full)
133 CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"
134 CODE+=createCode(f_interpolation, x, Q_r, loopindex=idx_full)
135 CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"
136 for i in xrange(len(Q_r_faces)):
137 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)
139 CODE+="\n} /* end of face %s */\n"%i
140 CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"
141 for i in xrange(len(Q_r_faces)):
142 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)
144 CODE+="\n} /* end of face %s */\n"%i
145 CODE+="\n} /* end of out_data_type branching */\n"
146 insertCode("Assemble_Interpolation_%sD.c"%DIM, { "SNIP" : CODE})
147
148 # gradient to Quad points
149 CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"
150 CODE+=createGradientCode(f_interpolation, x, Q, loopindex=idx_full)
151 CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"
152 CODE+=createGradientCode(f_interpolation, x, Q_r, loopindex=idx_full)
153 CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"
154 for i in xrange(len(Q_r_faces)):
155 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)
157 CODE+="\n} /* end of face %s */\n"%i
158 CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"
159 for i in xrange(len(Q_r_faces)):
160 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)
162 CODE+="\n} /* end of face %s */\n"%i
163 CODE+="\n} /* end of out_data_type branching */\n"
164 insertCode("Assemble_Gradient_%sD.c"%DIM, { "SNIP" : CODE})
165
166 #generate PDE assemblage
167 CODE,PRECODE = makePDE(S, x, Q, W, DIM=DIM, system=True)
168 insertCode("Assemble_PDE_System_%sD.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE } )
169 CODE,PRECODE = makePDE(S, x, Q_r, W_r, DIM=DIM, system=True)
170 insertCode("Assemble_PDE_System_%sD_reduced.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE })
171 CODE,PRECODE = makePDE(S, x, Q, W, DIM=DIM, system=False)
172 insertCode("Assemble_PDE_Single_%sD.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE })
173 CODE,PRECODE = makePDE(S, x, Q_r, W_r, DIM=DIM, system=False)
174 insertCode("Assemble_PDE_Single_%sD_reduced.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE })
175
176
177
178
179 def insertCode(fn, replacement):
180 TOP_LINE_FMT="/* GENERATOR %s TOP"
181 BOTTOM_LINE_FMT="/* GENERATOR %s BOTTOM"
182 TAG=" "
183
184 # read code and create back up
185 bkp_stream=open(fn+".bkp",'w')
186 old_code=""
187 for l in open(fn, 'r').readlines():
188 s=l.strip()
189 if len(s)>0:
190 bkp_stream.write(s+"\n")
191 old_code+=s+"\n"
192 bkp_stream.close()
193
194 for k in replacement.keys():
195 TOP_LINE=TOP_LINE_FMT%k
196 BOTTOM_LINE=BOTTOM_LINE_FMT%k
197 new_code=""
198 ignore=False
199 for l in old_code.split("\n"):
200 s=l.strip()
201 if len(s)>0:
202 if s.startswith(TOP_LINE):
203 ignore=True
204 new_code+=s+"\n"+replacement[k]
205 elif s.startswith(BOTTOM_LINE):
206 ignore=False
207 new_code+=s+"\n"
208 elif not ignore:
209 new_code+=s+"\n"
210 old_code=new_code
211 #insert tabs:
212 N=0
213 new_code=""
214 for l in old_code.split("\n"):
215 s=l.strip()
216 if s.find("/*")>-1 and s.find("*/")>-1 :
217 s2=(s[:s.find("/*")] + s[s.find("*/")+2:]).strip()
218 else:
219 s2=s
220 if len(s)>0:
221 if s2.startswith("}"):
222 N=max(N-1, 0)
223 new_code+=TAG*N + s + "\n"
224 if s2.endswith("{"):
225 N+=1
226 open(fn,'w').write(new_code)
227 print "file %s generated."%fn
228
229
230
231 def optimizeEvaluate(F, x, Q):
232 F2=[]
233 for f in F:
234 f0=[]
235 for q in xrange(len(Q)):
236 tmp=f
237 for i in range(len(x)): tmp=tmp.subs(x[i], Q[q][i])
238 f0.append(tmp)
239 F2.append(f0)
240
241 # collect arguments:
242 args=[]
243 for f in F2:
244 for fq in f:
245 for s in fq.atoms():
246 if isinstance(s, Symbol):
247 if collect(fq, s, evaluate=False).has_key(s):
248 if s.name.startswith("f_"):
249 if args.count(s) == 0 and not (collect(fq, s, evaluate=False)[s]).is_zero: args.append(s)
250 c0={}
251 cc=0
252 F_new=[]
253 for f in F2:
254 fq_new=[]
255 for fq in f:
256 s1={}
257 for a in args:
258 s=collect(fq, a, evaluate=False)
259 if s.has_key(a):
260 k0=None
261 for k, v in c0.items():
262 if s[a] == v: k0=k
263 if k0==None:
264 k0=Symbol("tmp0_%s"%cc)
265 c0[k0]=s[a]
266 cc+=1
267 if not s1.has_key(k0): s1[k0]=0
268 s1[k0] = s1[k0] + a
269 tmp=0
270 for k,v in s1.items():
271 tmp = tmp + k * v
272 fq_new.append(tmp)
273 F_new.append(fq_new)
274
275 return F_new, args, c0
276
277 def createGradientCode(F, x, Q, gridoffset="", loopindex=(0,1,2)):
278 DIM=len(loopindex)
279 dF=[ expand(diff(F, x[i])) for i in range(DIM)]
280 dF, args, consts = optimizeEvaluate(dF, x, Q)
281 k=""
282 N=""
283 M1=""
284 for i in xrange(len(Q[0])):
285 if len(k)==0:
286 k+="k%s"%i
287 N+="N%s"%i
288 M1+="M%s"%i
289 else:
290 k+=",k%s"%i
291 if i < DIM-1: N+=",N%s"%i
292 if i < DIM-1: M1+=",M%s"%i
293
294 k3=""
295 for i in xrange(DIM-1,-1,-1):
296 if loopindex[i] > -1: k3+=",k%i"%loopindex[i]
297 TXT=""
298 for a,v in consts.items():
299 TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))
300 TXT+="#pragma omp parallel for private(i%s)"%k3
301
302 for i in xrange(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])
304 TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"
305 for s in args:
306 k1=""
307 for i in xrange(DIM):
308 if s.name[2+i]=="0":
309 if loopindex[i]>-1:
310 k1+=",k%s"%i
311 elif loopindex[i] == -1:
312 k1+=",0"
313 elif loopindex[i] == -2:
314 k1+=",M%s-2"%i
315 else:
316 if loopindex[i]>-1:
317 k1+=",k%s+1"%i
318 elif loopindex[i] == -1:
319 k1+=",1"
320 elif loopindex[i] == -2:
321 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)
323 # interpolation to quadrature points
324 for q in xrange(len(Q)):
325 IDX2="INDEX%s(%s,%s)"%(DIM,k,N)
326 if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2
327 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]))
329 TXT+="} /* close component loop i */\n"
330 for i in xrange(DIM):
331 if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]
332 return TXT
333
334 def createCode(F, x, Q, gridoffset="", loopindex=(0,1,2)):
335 DIM=len(loopindex)
336 F, args, consts = optimizeEvaluate([F, ], x, Q)
337 F=F[0]
338 k=""
339 N=""
340 M1=""
341 for i in xrange(len(Q[0])):
342 if len(k)==0:
343 k+="k%s"%i
344 N+="N%s"%i
345 M1+="M%s"%i
346 else:
347 k+=",k%s"%i
348 if i < DIM-1: N+=",N%s"%i
349 if i < DIM-1: M1+=",M%s"%i
350
351 k3=""
352 for i in xrange(DIM-1,-1,-1):
353 if loopindex[i] > -1: k3+=",k%i"%loopindex[i]
354 TXT=""
355 for a,v in consts.items():
356 TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))
357 TXT+="#pragma omp parallel for private(i%s)"%k3
358 for i in xrange(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])
360 TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"
361 for s in args:
362 k1=""
363 for i in xrange(DIM):
364 if s.name[2+i]=="0":
365 if loopindex[i]>-1:
366 k1+=",k%s"%i
367 elif loopindex[i] == -1:
368 k1+=",0"
369 elif loopindex[i] == -2:
370 k1+=",M%s-2"%i
371 else:
372 if loopindex[i]>-1:
373 k1+=",k%s+1"%i
374 elif loopindex[i] == -1:
375 k1+=",1"
376 elif loopindex[i] == -2:
377 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)
379 # interpolation to quadrature points
380 for q in xrange(len(Q)):
381 IDX2="INDEX%s(%s,%s)"%(DIM,k,N)
382 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]))
384 TXT+="} /* close component loop i */\n"
385 for i in xrange(DIM):
386 if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]
387 return TXT
388
389 def subPoint(g,x, p):
390 out=g
391 for i in range(len(x)): out=out.subs(x[i], p[i])
392 # return out.evalf()
393 return out
394
395 def generatePDECode(DATA_A, EM,GLOBAL_TMP, system=False):
396 LOCAL_TMP={}
397 LOCAL2_TMP={}
398 OUT=""
399 for p in EM.keys():
400 # first we collect the constant coefficients in the expression (outside the parallel region):
401 tmp={}
402 for a in DATA_A:
403 c= collect(EM[p], a, evaluate=False)
404 if c.has_key(a):
405 if not GLOBAL_TMP.has_key(c[a]):
406 GLOBAL_TMP[c[a]]=Symbol("w%s"%(len(GLOBAL_TMP),))
407 if tmp.has_key(GLOBAL_TMP[c[a]]):
408 tmp[GLOBAL_TMP[c[a]]]+=a
409 else:
410 tmp[GLOBAL_TMP[c[a]]]=a
411 # now we find temporary values which are combinations of the input data:
412 em_new=[]
413 for t in tmp.keys():
414 if tmp[t] in DATA_A:
415 tt = t * tmp[t]
416 else:
417 if not LOCAL_TMP.has_key(tmp[t]):
418 LOCAL_TMP[tmp[t]]=Symbol("tmp%s_0"%len(LOCAL_TMP))
419 tt= t * LOCAL_TMP[tmp[t]]
420 if not LOCAL2_TMP.has_key(tt):
421 LOCAL2_TMP[tt]=Symbol("tmp%s_1"%len(LOCAL2_TMP))
422 em_new.append(LOCAL2_TMP[tt])
423 EM[p]=em_new
424
425 for p in EM.keys():
426 sum_em=0
427 for t in EM[p]: sum_em=sum_em + t
428 if isinstance(p, int) :
429 if system:
430 OUT+=" EM_F[INDEX2(k,%s,p.numEqu)]+=%s;\n"%(p,ccode(sum_em))
431 else:
432 OUT+=" EM_F[%s]+=%s;\n"%(p,ccode(sum_em))
433 else:
434 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))
436 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))
438
439 OUT2=""
440 for p in LOCAL_TMP:
441 OUT2+=" const register double %s = %s;\n"%(LOCAL_TMP[p],ccode(p))
442 for p in LOCAL2_TMP:
443 OUT2+=" const register double %s = %s;\n"%(LOCAL2_TMP[p],ccode(p))
444 return OUT2+OUT
445
446 def makePDE(S, x, Q, W, DIM=2, system=False):
447 GLOBAL_TMP={}
448 GLOBAL_N=0
449 PRECODE=""
450 CODE="""
451 /**************************************************************/
452 /* process A: */
453 /**************************************************************/
454 if (NULL!=A_p) {
455 add_EM_S=TRUE;
456 """
457 if len(Q) > 1:
458 CODE+="if (extendedA) {\n"
459 if system: CODE+= """for (k=0;k<p.numEqu;k++) {
460 for (m=0;m<p.numComp;m++) {
461 """
462 DATA_A=[]
463 CODE2=""
464 EM = {}
465 for i in range(len(S)):
466 for j in range(len(S)):
467 EM[(i,j)] = 0
468
469 for q in xrange(len(Q)):
470 for di in range(DIM):
471 for dj in range(DIM):
472 A_name="A_%d%d_%d"%(di,dj,q)
473 A=Symbol(A_name)
474 DATA_A.append(A)
475 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)
477 else:
478 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)):
480 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) ] )
482 CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system)
483
484 if system: CODE+="}\n }\n"
485 CODE+="} else { /* constant data */\n"
486 if system:
487 CODE+= """for (k=0;k<p.numEqu;k++) {
488 for (m=0;m<p.numComp;m++) {
489 """
490 DATA_A=[]
491 CODE2=""
492 EM = {}
493 for i in range(len(S)):
494 for j in range(len(S)):
495 EM[(i,j)] = 0
496 for di in range(DIM):
497 for dj in range(DIM):
498 A_name="A_%d%d"%(di,dj)
499 A=Symbol(A_name)
500 DATA_A.append(A)
501 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)
503 else:
504 CODE2+="const register double %s = A_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, dj, DIM)
505 for q in xrange(len(Q)):
506 for i in range(len(S)):
507 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) ] )
509 CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system)
510 if system: CODE+="}\n }\n"
511 if len(Q) > 1: CODE+="}\n"
512 CODE+="}\n"
513 #BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
514 CODE+="""
515 /**************************************************************/
516 /* process B: */
517 /**************************************************************/
518 if (NULL!=B_p) {
519 add_EM_S=TRUE;
520 """
521 if len(Q) > 1:
522 CODE+="if (extendedB) {\n"
523 if system: CODE+= """for (k=0;k<p.numEqu;k++) {
524 for (m=0;m<p.numComp;m++) {
525 """
526 DATA_B=[]
527 CODE2=""
528 EM = {}
529 for i in range(len(S)):
530 for j in range(len(S)):
531 EM[(i,j)] = 0
532
533 for q in xrange(len(Q)):
534 for di in range(DIM):
535 A_name="B_%d_%d"%(di,q)
536 A=Symbol(A_name)
537 DATA_B.append(A)
538 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)
540 else:
541 CODE2+="const register double %s = B_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, q, DIM)
542 for i in range(len(S)):
543 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) ] )
545 CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system)
546
547 if system: CODE+="}\n }\n"
548 CODE+="} else { /* constant data */\n"
549 if system:
550 CODE+= """for (k=0;k<p.numEqu;k++) {
551 for (m=0;m<p.numComp;m++) {
552 """
553 DATA_B=[]
554 CODE2=""
555 EM = {}
556 for i in range(len(S)):
557 for j in range(len(S)):
558 EM[(i,j)] = 0
559 for di in range(DIM):
560 A_name="B_%d"%(di)
561 A=Symbol(A_name)
562 DATA_B.append(A)
563 if system:
564 CODE2+="const register double %s = B_p[INDEX3(k,%s,m, p.numEqu,%s)];\n"%(A_name, di, DIM)
565 else:
566 CODE2+="const register double %s = B_p[%s];\n"%(A_name, di)
567 for q in xrange(len(Q)):
568 for i in range(len(S)):
569 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) ] )
571 CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system)
572 if system: CODE+="}\n }\n"
573 if len(Q) > 1: CODE+="}\n"
574 CODE+="}\n"
575 #CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
576 CODE+="""
577 /**************************************************************/
578 /* process C: */
579 /**************************************************************/
580 if (NULL!=C_p) {
581 add_EM_S=TRUE;
582 """
583 if len(Q) > 1:
584 CODE+="if (extendedC) {\n"
585 if system: CODE+= """for (k=0;k<p.numEqu;k++) {
586 for (m=0;m<p.numComp;m++) {
587 """
588 DATA_C=[]
589 CODE2=""
590 EM = {}
591 for i in range(len(S)):
592 for j in range(len(S)):
593 EM[(i,j)] = 0
594
595 for q in xrange(len(Q)):
596 for dj in range(DIM):
597 A_name="C_%d_%d"%(dj,q)
598 A=Symbol(A_name)
599 DATA_C.append(A)
600 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)
602 else:
603 CODE2+="const register double %s = C_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj, q, DIM)
604 for i in range(len(S)):
605 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) ] )
607 CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system)
608
609 if system: CODE+="}\n }\n"
610 CODE+="} else { /* constant data */\n"
611 if system:
612 CODE+= """for (k=0;k<p.numEqu;k++) {
613 for (m=0;m<p.numComp;m++) {
614 """
615 DATA_C=[]
616 CODE2=""
617 EM = {}
618 for i in range(len(S)):
619 for j in range(len(S)):
620 EM[(i,j)] = 0
621 for dj in range(DIM):
622 A_name="C_%d"%(dj)
623 A=Symbol(A_name)
624 DATA_C.append(A)
625 if system:
626 CODE2+="const register double %s = C_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, dj)
627 else:
628 CODE2+="const register double %s = C_p[%s];\n"%(A_name, dj)
629 for q in xrange(len(Q)):
630 for i in range(len(S)):
631 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) ] )
633 CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system)
634 if system: CODE+="}\n }\n"
635 if len(Q) > 1: CODE+="}\n"
636 CODE+="}\n"
637 #DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD
638 CODE+="""
639 /**************************************************************/
640 /* process D: */
641 /**************************************************************/
642 if (NULL!=D_p) {
643 add_EM_S=TRUE;
644 """
645 if len(Q) > 1:
646 CODE+="if (extendedD) {\n"
647 if system: CODE+= """for (k=0;k<p.numEqu;k++) {
648 for (m=0;m<p.numComp;m++) {
649 """
650 DATA_D=[]
651 CODE2=""
652 EM = {}
653 for i in range(len(S)):
654 for j in range(len(S)):
655 EM[(i,j)] = 0
656
657 for q in xrange(len(Q)):
658 A_name="D_%d"%(q,)
659 A=Symbol(A_name)
660 DATA_D.append(A)
661 if system:
662 CODE2+="const register double %s = D_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, q)
663 else:
664 CODE2+="const register double %s = D_p[%s];\n"%(A_name, q)
665 for i in range(len(S)):
666 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) ] )
668 CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system)
669
670 if system: CODE+="}\n }\n"
671 CODE+="} else { /* constant data */\n"
672 if system:
673 CODE+= """for (k=0;k<p.numEqu;k++) {
674 for (m=0;m<p.numComp;m++) {
675 """
676 DATA_D=[]
677 CODE2=""
678 EM = {}
679 for i in range(len(S)):
680 for j in range(len(S)):
681 EM[(i,j)] = 0
682 if 1:
683 A_name="D_0"
684 A=Symbol(A_name)
685 DATA_D.append(A)
686 if system:
687 CODE2+="const register double %s = D_p[INDEX2(k,m, p.numEqu)];\n"%(A_name,)
688 else:
689 CODE2+="const register double %s = D_p[0];\n"%(A_name,)
690 for q in xrange(len(Q)):
691 for i in range(len(S)):
692 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) ] )
694 CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system)
695 if system: CODE+="}\n }\n"
696 if len(Q) > 1: CODE+="}\n"
697 CODE+="}\n"
698
699 #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
700 CODE+="""
701 /**************************************************************/
702 /* process X: */
703 /**************************************************************/
704 if (NULL!=X_p) {
705 add_EM_F=TRUE;
706 """
707 if len(Q) > 1:
708 CODE+="if (extendedX) {\n"
709 if system: CODE+= """for (k=0;k<p.numEqu;k++) {
710 """
711 DATA_X=[]
712 CODE2=""
713 EM = {}
714 for i in range(len(S)):
715 EM[i] = 0
716
717 for q in xrange(len(Q)):
718 for dj in range(DIM):
719 A_name="X_%d_%d"%(dj,q)
720 A=Symbol(A_name)
721 DATA_X.append(A)
722 if system:
723 CODE2+="const register double %s = X_p[INDEX3(k,%s, %s, p.numEqu,%s)];\n"%(A_name, dj, q, DIM)
724 else:
725 CODE2+="const register double %s = X_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj,q,DIM)
726 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) ] )
728 CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system)
729
730 if system: CODE+="}\n"
731 CODE+="} else { /* constant data */\n"
732 if system:
733 CODE+= """for (k=0;k<p.numEqu;k++) {
734 """
735 DATA_X=[]
736 CODE2=""
737 EM = {}
738 for i in range(len(S)):
739 EM[i] = 0
740 for dj in range(DIM):
741 A_name="X_%d"%(dj)
742 A=Symbol(A_name)
743 DATA_X.append(A)
744 if system:
745 CODE2+="const register double %s = X_p[INDEX2(k,%s, p.numEqu)];\n"%(A_name, dj)
746 else:
747 CODE2+="const register double %s = X_p[%s];\n"%(A_name, dj)
748 for q in xrange(len(Q)):
749 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) ] )
751 CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system)
752 if system: CODE+="}\n"
753 if len(Q) > 1: CODE+="}\n"
754 CODE+="}\n"
755
756 #YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYy
757 CODE+="""
758 /**************************************************************/
759 /* process Y: */
760 /**************************************************************/
761 if (NULL!=Y_p) {
762 add_EM_F=TRUE;
763 """
764 if len(Q) > 1:
765 CODE+="if (extendedY) {\n"
766 if system: CODE+= """for (k=0;k<p.numEqu;k++) {
767 """
768 DATA_Y=[]
769 CODE2=""
770 EM = {}
771 for i in range(len(S)):
772 EM[i] = 0
773
774 for q in xrange(len(Q)):
775 A_name="Y_%d"%(q,)
776 A=Symbol(A_name)
777 DATA_Y.append(A)
778 if system:
779 CODE2+="const register double %s = Y_p[INDEX3(k,%s, p.numEqu)];\n"%(A_name, q)
780 else:
781 CODE2+="const register double %s = Y_p[%s];\n"%(A_name, q)
782 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) ] )
784 CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system)
785
786 if system: CODE+="}\n"
787 CODE+="} else { /* constant data */\n"
788 if system:
789 CODE+= """for (k=0;k<p.numEqu;k++) {
790 """
791 DATA_Y=[]
792 CODE2=""
793 EM = {}
794 for i in range(len(S)):
795 EM[i] = 0
796 if 1:
797 A_name="Y_0"
798 A=Symbol(A_name)
799 DATA_Y.append(A)
800 if system:
801 CODE2+="const register double %s = Y_p[k];\n"%(A_name,)
802 else:
803 CODE2+="const register double %s = Y_p[0];\n"%(A_name,)
804 for q in xrange(len(Q)):
805 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) ] )
807 CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system)
808 if system: CODE+="}\n"
809 if len(Q) > 1: CODE+="}\n"
810 CODE+="}\n"
811
812 for t in GLOBAL_TMP:
813 PRECODE+="const double %s = %s;\n"%(GLOBAL_TMP[t],ccode(t.evalf(n=DIGITS)))
814 return CODE, PRECODE
815
816 for d in [3]:
817 generate(d)

  ViewVC Help
Powered by ViewVC 1.1.26