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

  ViewVC Help
Powered by ViewVC 1.1.26