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

  ViewVC Help
Powered by ViewVC 1.1.26