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

  ViewVC Help
Powered by ViewVC 1.1.26