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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 gross 3543 """
2 caltinay 3645 this script generates the assemblage routine for the ripley rectangular grid solver
3 gross 3543
4    
5     """
6     from sympy import *
7    
8     DIGITS=20
9    
10     #
11     # quadrature form in [0,1]
12     #
13 gross 3689
14     q0=1/sqrt(RealNumber(3))
15     qD1=[ (1-q0)/2, (1+q0)/2 ]
16     wD1=[RealNumber(1)/2, RealNumber(1)/2 ]
17    
18 caltinay 3792
19     q0=1/sqrt(RealNumber(3))
20     qD1=[ (1-q0)/2, (1+q0)/2 ]
21     wD1=[RealNumber(1)/2, RealNumber(1)/2 ]
22    
23 gross 3543 qD1_r=[ RealNumber(1)/2 ]
24     wD1_r= [RealNumber(1)]
25 gross 3689
26 gross 3543 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 caltinay 3635 CODE+="\n} /* end of out_data_type branching */\n"
146 gross 3641 insertCode("Assemble_Interpolation_%sD.c"%DIM, { "SNIP" : CODE})
147 caltinay 3792
148     # gradient to Quad points
149 gross 3543 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 caltinay 3635 CODE+="\n} /* end of out_data_type branching */\n"
164 gross 3641 insertCode("Assemble_Gradient_%sD.c"%DIM, { "SNIP" : CODE})
165 caltinay 3792
166 gross 3641 #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 gross 3543
178    
179 gross 3641 def insertCode(fn, replacement):
180     TOP_LINE_FMT="/* GENERATOR %s TOP"
181     BOTTOM_LINE_FMT="/* GENERATOR %s BOTTOM"
182 gross 3543 TAG=" "
183    
184 gross 3641 # read code and create back up
185 gross 3543 bkp_stream=open(fn+".bkp",'w')
186 gross 3641 old_code=""
187 gross 3543 for l in open(fn, 'r').readlines():
188     s=l.strip()
189     if len(s)>0:
190     bkp_stream.write(s+"\n")
191 gross 3641 old_code+=s+"\n"
192 gross 3543 bkp_stream.close()
193 gross 3641
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 gross 3543 N=0
213 gross 3641 new_code=""
214     for l in old_code.split("\n"):
215 gross 3543 s=l.strip()
216 gross 3641 if s.find("/*")>-1 and s.find("*/")>-1 :
217     s2=(s[:s.find("/*")] + s[s.find("*/")+2:]).strip()
218     else:
219     s2=s
220 gross 3543 if len(s)>0:
221 gross 3641 if s2.startswith("}"):
222 gross 3543 N=max(N-1, 0)
223 gross 3641 new_code+=TAG*N + s + "\n"
224     if s2.endswith("{"):
225 gross 3543 N+=1
226 gross 3641 open(fn,'w').write(new_code)
227     print "file %s generated."%fn
228 gross 3543
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 gross 3641 TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)
323 gross 3543 # 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 gross 3641 TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)
379 gross 3543 # 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 gross 3641
389 gross 3543 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 gross 3641 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 gross 3543
439 gross 3641 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 caltinay 3645 CODE2+="const register double %s = X_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj,q,DIM)
726 gross 3641 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 caltinay 3792 for d in [3]:
817 gross 3543 generate(d)

  ViewVC Help
Powered by ViewVC 1.1.26