/[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 3635 - (hide annotations)
Fri Oct 21 00:26:01 2011 UTC (8 years, 1 month ago) by caltinay
File MIME type: text/x-python
File size: 12589 byte(s)
Fixed minor issues.

1 gross 3543 """
2     this script generates the assamblage 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     qD1_r=[ RealNumber(1)/2 ]
17     wD1_r= [RealNumber(1)]
18     print "1D quadrature nodes =",qD1
19     print "1D quadrature weights =",wD1
20     print "1D reduced quadrature nodes =",qD1_r
21     print "1D reduced quadrature weights =",wD1_r
22    
23     def generate(DIM):
24     # variables
25     h=[Symbol('h%s'%i) for i in xrange(DIM) ]
26     x=[Symbol('x%s'%i) for i in xrange(DIM) ]
27     # shape functions: labeling differently from finley!!!!
28     S=[ 1 ]
29     GridOffset=[[]]
30     GridOffsetText=[ "" ]
31     idx_full=[]
32     Q = [ [] ]
33     W = [ 1 ]
34     Q_r = [ [] ]
35     W_r = [ 1 ]
36     for i in xrange(DIM):
37     idx_full.append(i)
38     S = [ s * (1-x[i]/h[i]) for s in S] + [ s * x[i]/h[i] for s in S]
39     GridOffset = [ s + [0] for s in GridOffset] + [ s + [1] for s in GridOffset]
40     GridOffsetText = [ s + "0" for s in GridOffsetText] + [ s + "1" for s in GridOffsetText]
41    
42     # generate quadrature points in element
43     Qnew, Wnew =[], []
44     for j in xrange(len(qD1)):
45     Qnew += [ q + [qD1[j]*h[i],] for q in Q ]
46     Wnew += [ w * wD1[j]*h[i] for w in W ]
47     Q, W=Qnew, Wnew
48    
49     Qnew, Wnew =[], []
50     for j in xrange(len(qD1_r)):
51     Qnew += [ q + [qD1_r[j]*h[i],] for q in Q_r ]
52     Wnew += [ w * wD1_r[j]*h[i] for w in W_r ]
53     Q_r, W_r=Qnew, Wnew
54    
55     # now the same thing for the faces:
56     Q_faces=[]
57     W_faces=[]
58     Q_r_faces=[]
59     W_r_faces=[]
60     idx_faces=[]
61     for face in xrange(DIM):
62     Q_left, W_left = [ [] ], [ 1 ]
63     Q_left_r, W_left_r = [ [] ], [ 1 ]
64     Q_right, W_right = [ [] ], [ 1 ]
65     Q_right_r, W_right_r = [ [] ], [ 1 ]
66     idx_left,idx_right=[],[]
67     for i in xrange(DIM):
68     # generate quadrature points in element
69     Q_left_new, W_left_new =[], []
70     Q_right_new, W_right_new =[], []
71    
72     if face == i :
73     idx_left.append(-1)
74     idx_right.append(-2)
75     Q_left_new += [ q + [0.,] for q in Q_left ]
76     W_left_new += [ w * 1. for w in W_left ]
77     Q_right_new += [ q + [ h[i], ] for q in Q_right ]
78     W_right_new += [ w * 1. for w in W_right ]
79     else:
80     idx_left.append(i)
81     idx_right.append(i)
82     for j in xrange(len(qD1)):
83     Q_left_new += [ q + [qD1[j]*h[i],] for q in Q_left ]
84     W_left_new += [ w * wD1[j]*h[i] for w in W_left ]
85     Q_right_new += [ q + [qD1[j]*h[i],] for q in Q_right ]
86     W_right_new += [ w * wD1[j]*h[i] for w in W_right ]
87     Q_left, W_left=Q_left_new, W_left_new
88     Q_right, W_right=Q_right_new, W_right_new
89     Q_faces=Q_faces+[Q_left, Q_right]
90     W_faces=W_faces+[W_left, W_right]
91     idx_faces=idx_faces + [ idx_left, idx_right ]
92    
93     Q_left, W_left = [ [] ], [ 1 ]
94     Q_left_r, W_left_r = [ [] ], [ 1 ]
95     Q_right, W_right = [ [] ], [ 1 ]
96     Q_right_r, W_right_r = [ [] ], [ 1 ]
97     for i in xrange(DIM):
98     # generate quadrature points in element
99     Q_left_new, W_left_new =[], []
100     Q_right_new, W_right_new =[], []
101     if face == i :
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     for j in xrange(len(qD1_r)):
108     Q_left_new += [ q + [qD1_r[j]*h[i],] for q in Q_left ]
109     W_left_new += [ w * wD1_r[j]*h[i] for w in W_left ]
110     Q_right_new += [ q + [qD1_r[j]*h[i],] for q in Q_right ]
111     W_right_new += [ w * wD1_r[j]*h[i] for w in W_right ]
112     Q_left, W_left=Q_left_new, W_left_new
113     Q_right, W_right=Q_right_new, W_right_new
114     Q_r_faces=Q_r_faces+[Q_left, Q_right]
115     W_r_faces=W_r_faces+[W_left, W_right]
116    
117     # the interpolation function
118     f_interpolation = 0
119     for i in xrange(len(GridOffsetText)):
120     f_interpolation = f_interpolation + S[i] * Symbol("f_%s"%GridOffsetText[i])
121    
122     # interpolation to Quad points
123     CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"
124     CODE+=createCode(f_interpolation, x, Q, loopindex=idx_full)
125     CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"
126     CODE+=createCode(f_interpolation, x, Q_r, loopindex=idx_full)
127     CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"
128     for i in xrange(len(Q_r_faces)):
129     CODE+="if (face_offset(%s)>-1) {\n"%i
130     CODE+=createCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i)
131     CODE+="\n} /* end of face %s */\n"%i
132     CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"
133     for i in xrange(len(Q_r_faces)):
134     CODE+="if (face_offset(%s)>-1) {\n"%i
135     CODE+=createCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i)
136     CODE+="\n} /* end of face %s */\n"%i
137 caltinay 3635 CODE+="\n} /* end of out_data_type branching */\n"
138 gross 3543 insertCode("Assemble_Interpolation_%sD.c"%DIM, CODE)
139    
140    
141     # interpolation to Quad points
142     CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"
143     CODE+=createGradientCode(f_interpolation, x, Q, loopindex=idx_full)
144     CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"
145     CODE+=createGradientCode(f_interpolation, x, Q_r, loopindex=idx_full)
146     CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"
147     for i in xrange(len(Q_r_faces)):
148     CODE+="if (face_offset(%s)>-1) {\n"%i
149     CODE+=createGradientCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i)
150     CODE+="\n} /* end of face %s */\n"%i
151     CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"
152     for i in xrange(len(Q_r_faces)):
153     CODE+="if (face_offset(%s)>-1) {\n"%i
154     CODE+=createGradientCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i)
155     CODE+="\n} /* end of face %s */\n"%i
156 caltinay 3635 CODE+="\n} /* end of out_data_type branching */\n"
157 gross 3543 insertCode("Assemble_Gradient_%sD.c"%DIM, CODE)
158    
159    
160     def insertCode(fn, code):
161     TOP_LINE="/* GENERATOR SNIP TOP"
162     BOTTOM_LINE="/* GENERATOR SNIP BOTTOM"
163     TAG=" "
164    
165     bkp_stream=open(fn+".bkp",'w')
166     new_code=""
167     ignore=False
168     for l in open(fn, 'r').readlines():
169     s=l.strip()
170     if len(s)>0:
171     bkp_stream.write(s+"\n")
172     if s.startswith(TOP_LINE):
173     ignore=True
174     new_code+=s+"\n"+code
175     elif s.startswith(BOTTOM_LINE):
176     ignore=False
177     new_code+=s+"\n"
178     elif not ignore:
179     new_code+=s+"\n"
180     bkp_stream.close()
181     N=0
182     new_code2=""
183     for l in new_code.split("\n"):
184     s=l.strip()
185     if len(s)>0:
186     if s.startswith("}"):
187     N=max(N-1, 0)
188     new_code2+=TAG*N + s + "\n"
189     if s.endswith("{"):
190     N+=1
191     open(fn,'w').write(new_code2)
192    
193    
194    
195     def optimizeEvaluate(F, x, Q):
196     F2=[]
197     for f in F:
198     f0=[]
199     for q in xrange(len(Q)):
200     tmp=f
201     for i in range(len(x)): tmp=tmp.subs(x[i], Q[q][i])
202     f0.append(tmp)
203     F2.append(f0)
204    
205     # collect arguments:
206     args=[]
207     for f in F2:
208     for fq in f:
209     for s in fq.atoms():
210     if isinstance(s, Symbol):
211     if collect(fq, s, evaluate=False).has_key(s):
212     if s.name.startswith("f_"):
213     if args.count(s) == 0 and not (collect(fq, s, evaluate=False)[s]).is_zero: args.append(s)
214     c0={}
215     cc=0
216     F_new=[]
217     for f in F2:
218     fq_new=[]
219     for fq in f:
220     s1={}
221     for a in args:
222     s=collect(fq, a, evaluate=False)
223     if s.has_key(a):
224     k0=None
225     for k, v in c0.items():
226     if s[a] == v: k0=k
227     if k0==None:
228     k0=Symbol("tmp0_%s"%cc)
229     c0[k0]=s[a]
230     cc+=1
231     if not s1.has_key(k0): s1[k0]=0
232     s1[k0] = s1[k0] + a
233     tmp=0
234     for k,v in s1.items():
235     tmp = tmp + k * v
236     fq_new.append(tmp)
237     F_new.append(fq_new)
238    
239     return F_new, args, c0
240    
241     def createGradientCode(F, x, Q, gridoffset="", loopindex=(0,1,2)):
242     DIM=len(loopindex)
243     dF=[ expand(diff(F, x[i])) for i in range(DIM)]
244     dF, args, consts = optimizeEvaluate(dF, x, Q)
245     k=""
246     N=""
247     M1=""
248     for i in xrange(len(Q[0])):
249     if len(k)==0:
250     k+="k%s"%i
251     N+="N%s"%i
252     M1+="M%s"%i
253     else:
254     k+=",k%s"%i
255     if i < DIM-1: N+=",N%s"%i
256     if i < DIM-1: M1+=",M%s"%i
257    
258     k3=""
259     for i in xrange(DIM-1,-1,-1):
260     if loopindex[i] > -1: k3+=",k%i"%loopindex[i]
261     TXT=""
262     for a,v in consts.items():
263     TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))
264     TXT+="#pragma omp parallel for private(i%s)"%k3
265    
266     for i in xrange(DIM-1,-1,-1):
267     if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i])
268     TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"
269     for s in args:
270     k1=""
271     for i in xrange(DIM):
272     if s.name[2+i]=="0":
273     if loopindex[i]>-1:
274     k1+=",k%s"%i
275     elif loopindex[i] == -1:
276     k1+=",0"
277     elif loopindex[i] == -2:
278     k1+=",M%s-2"%i
279     else:
280     if loopindex[i]>-1:
281     k1+=",k%s+1"%i
282     elif loopindex[i] == -1:
283     k1+=",1"
284     elif loopindex[i] == -2:
285     k1+=",M%s-1"%i
286     TXT+="register const double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)
287     # interpolation to quadrature points
288     for q in xrange(len(Q)):
289     IDX2="INDEX%s(%s,%s)"%(DIM,k,N)
290     if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2
291     for i in range(DIM):
292     TXT+="out[INDEX4(i,%s,%s,%s,NCOMP,%s,%s)] = %s;\n"%(i,q,IDX2,DIM,len(Q),ccode(dF[i][q]))
293     TXT+="} /* close component loop i */\n"
294     for i in xrange(DIM):
295     if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]
296     return TXT
297    
298     def createCode(F, x, Q, gridoffset="", loopindex=(0,1,2)):
299     DIM=len(loopindex)
300     F, args, consts = optimizeEvaluate([F, ], x, Q)
301     F=F[0]
302     k=""
303     N=""
304     M1=""
305     for i in xrange(len(Q[0])):
306     if len(k)==0:
307     k+="k%s"%i
308     N+="N%s"%i
309     M1+="M%s"%i
310     else:
311     k+=",k%s"%i
312     if i < DIM-1: N+=",N%s"%i
313     if i < DIM-1: M1+=",M%s"%i
314    
315     k3=""
316     for i in xrange(DIM-1,-1,-1):
317     if loopindex[i] > -1: k3+=",k%i"%loopindex[i]
318     TXT=""
319     for a,v in consts.items():
320     TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))
321     TXT+="#pragma omp parallel for private(i%s)"%k3
322     for i in xrange(DIM-1,-1,-1):
323     if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i])
324     TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"
325     for s in args:
326     k1=""
327     for i in xrange(DIM):
328     if s.name[2+i]=="0":
329     if loopindex[i]>-1:
330     k1+=",k%s"%i
331     elif loopindex[i] == -1:
332     k1+=",0"
333     elif loopindex[i] == -2:
334     k1+=",M%s-2"%i
335     else:
336     if loopindex[i]>-1:
337     k1+=",k%s+1"%i
338     elif loopindex[i] == -1:
339     k1+=",1"
340     elif loopindex[i] == -2:
341     k1+=",M%s-1"%i
342     TXT+="register const double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)
343     # interpolation to quadrature points
344     for q in xrange(len(Q)):
345     IDX2="INDEX%s(%s,%s)"%(DIM,k,N)
346     if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2
347     TXT+="out[INDEX3(i,%s,%s,NCOMP,%s)] = %s;\n"%(q,IDX2,len(Q),ccode(F[q]))
348     TXT+="} /* close component loop i */\n"
349     for i in xrange(DIM):
350     if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]
351     return TXT
352    
353     def subPoint(g,x, p):
354     out=g
355     for i in range(len(x)): out=out.subs(x[i], p[i])
356     # return out.evalf()
357     return out
358    
359    
360     for d in [2,3 ]:
361     generate(d)

  ViewVC Help
Powered by ViewVC 1.1.26