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

Contents of /trunk/ripley/src/generate_assemblage_cpp.py

Parent Directory Parent Directory | Revision Log Revision Log


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

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=1/sqrt(RealNumber(3))
14 qD1=[ (1-q0)/2, (1+q0)/2 ]
15 wD1=[RealNumber(1)/2, RealNumber(1)/2 ]
16
17 qD1_r=[ RealNumber(1)/2 ]
18 wD1_r= [RealNumber(1)]
19
20 print "1D quadrature nodes =",qD1
21 print "1D quadrature weights =",wD1
22 print "1D reduced quadrature nodes =",qD1_r
23 print "1D reduced quadrature weights =",wD1_r
24
25 def generate(DIM, filename):
26 # variables
27 h=[Symbol('h%s'%i) for i in xrange(DIM) ]
28 x=[Symbol('x%s'%i) for i in xrange(DIM) ]
29 # shape functions: labeling differently from finley!!!!
30 S=[ 1 ]
31 GridOffset=[[]]
32 GridOffsetText=[ "" ]
33 idx_full=[]
34 Q = [ [] ]
35 W = [ 1 ]
36 Q_r = [ [] ]
37 W_r = [ 1 ]
38 for i in xrange(DIM):
39 idx_full.append(i)
40 S = [ s * (1-x[i]/h[i]) for s in S] + [ s * x[i]/h[i] for s in S]
41 GridOffset = [ s + [0] for s in GridOffset] + [ s + [1] for s in GridOffset]
42 GridOffsetText = [ s + "0" for s in GridOffsetText] + [ s + "1" for s in GridOffsetText]
43
44 # generate quadrature points in element
45 Qnew, Wnew =[], []
46 for j in xrange(len(qD1)):
47 Qnew += [ q + [qD1[j]*h[i],] for q in Q ]
48 Wnew += [ w * wD1[j]*h[i] for w in W ]
49 Q, W=Qnew, Wnew
50
51 Qnew, Wnew =[], []
52 for j in xrange(len(qD1_r)):
53 Qnew += [ q + [qD1_r[j]*h[i],] for q in Q_r ]
54 Wnew += [ w * wD1_r[j]*h[i] for w in W_r ]
55 Q_r, W_r=Qnew, Wnew
56
57 # now the same thing for the faces:
58 Q_faces=[]
59 W_faces=[]
60 Q_r_faces=[]
61 W_r_faces=[]
62 idx_faces=[]
63 for face in xrange(DIM):
64 Q_left, W_left = [ [] ], [ 1 ]
65 Q_left_r, W_left_r = [ [] ], [ 1 ]
66 Q_right, W_right = [ [] ], [ 1 ]
67 Q_right_r, W_right_r = [ [] ], [ 1 ]
68 idx_left,idx_right=[],[]
69 for i in xrange(DIM):
70 # generate quadrature points in element
71 Q_left_new, W_left_new =[], []
72 Q_right_new, W_right_new =[], []
73
74 if face == i :
75 idx_left.append(-1)
76 idx_right.append(-2)
77 Q_left_new += [ q + [0.,] for q in Q_left ]
78 W_left_new += [ w * 1. for w in W_left ]
79 Q_right_new += [ q + [ h[i], ] for q in Q_right ]
80 W_right_new += [ w * 1. for w in W_right ]
81 else:
82 idx_left.append(i)
83 idx_right.append(i)
84 for j in xrange(len(qD1)):
85 Q_left_new += [ q + [qD1[j]*h[i],] for q in Q_left ]
86 W_left_new += [ w * wD1[j]*h[i] for w in W_left ]
87 Q_right_new += [ q + [qD1[j]*h[i],] for q in Q_right ]
88 W_right_new += [ w * wD1[j]*h[i] for w in W_right ]
89 Q_left, W_left=Q_left_new, W_left_new
90 Q_right, W_right=Q_right_new, W_right_new
91 Q_faces=Q_faces+[Q_left, Q_right]
92 W_faces=W_faces+[W_left, W_right]
93 idx_faces=idx_faces + [ idx_left, idx_right ]
94
95 Q_left, W_left = [ [] ], [ 1 ]
96 Q_left_r, W_left_r = [ [] ], [ 1 ]
97 Q_right, W_right = [ [] ], [ 1 ]
98 Q_right_r, W_right_r = [ [] ], [ 1 ]
99 for i in xrange(DIM):
100 # generate quadrature points in element
101 Q_left_new, W_left_new =[], []
102 Q_right_new, W_right_new =[], []
103 if face == i :
104 Q_left_new += [ q + [0.,] for q in Q_left ]
105 W_left_new += [ w * 1. for w in W_left ]
106 Q_right_new += [ q + [ h[i], ] for q in Q_right ]
107 W_right_new += [ w * 1. for w in W_right ]
108 else:
109 for j in xrange(len(qD1_r)):
110 Q_left_new += [ q + [qD1_r[j]*h[i],] for q in Q_left ]
111 W_left_new += [ w * wD1_r[j]*h[i] for w in W_left ]
112 Q_right_new += [ q + [qD1_r[j]*h[i],] for q in Q_right ]
113 W_right_new += [ w * wD1_r[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_r_faces=Q_r_faces+[Q_left, Q_right]
117 W_r_faces=W_r_faces+[W_left, W_right]
118
119 # the interpolation function
120 f_interpolation = 0
121 for i in xrange(len(GridOffsetText)):
122 f_interpolation = f_interpolation + S[i] * Symbol("f_%s"%GridOffsetText[i])
123
124 # integration to Quad points
125 #CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"
126 #CODE+=createIntegrationCode(Q, W, loopindex=idx_full)
127 #CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"
128 #CODE+=createIntegrationCode(Q_r, W_r, loopindex=idx_full)
129 #CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"
130 #for i in xrange(len(Q_faces)):
131 # CODE+="if (face_offset(%s)>-1) {\n"%i
132 # CODE+=createIntegrationCode(Q_faces[i], W_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i)
133 # CODE+="\n} /* end of face %s */\n"%i
134 #CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"
135 #for i in xrange(len(Q_r_faces)):
136 # CODE+="if (face_offset(%s)>-1) {\n"%i
137 # CODE+=createIntegrationCode(Q_r_faces[i], W_r_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i)
138 # CODE+="\n} /* end of face %s */\n"%i
139 #CODE+="\n} /* end of out_data_type branching */\n"
140 #insertCode("Assemble_Integration_%sD.c"%DIM, { "SNIP" : CODE})
141 #1/0
142
143 # interpolation to Quad points
144 #CODE=createCode(f_interpolation, x, Q, loopindex=idx_full)
145 #insertCode(filename, { "SNIP_INTERPOLATE_ELEMENTS" : CODE})
146 #CODE=createCode(f_interpolation, x, Q_r, loopindex=idx_full)
147 #insertCode(filename, { "SNIP_INTERPOLATE_REDUCED_ELEMENTS" : CODE})
148 #CODE=""
149 #for i in xrange(len(Q_faces)):
150 # CODE+="if (m_faceOffset[%s] > -1) {\n"%i
151 # CODE+=createCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="m_faceOffset[%s]"%i)
152 # CODE+="\n} /* end of face %s */\n"%i
153 #insertCode(filename, { "SNIP_INTERPOLATE_FACES" : CODE})
154 #CODE=""
155 #for i in xrange(len(Q_r_faces)):
156 # CODE+="if (m_faceOffset[%s] > -1) {\n"%i
157 # CODE+=createCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="m_faceOffset[%s]"%i)
158 # CODE+="\n} /* end of face %s */\n"%i
159 #insertCode(filename, { "SNIP_INTERPOLATE_REDUCED_FACES" : CODE})
160
161 ## gradient to Quad points
162 #CODE=createGradientCode(f_interpolation, x, Q, loopindex=idx_full)
163 #insertCode(filename, { "SNIP_GRAD_ELEMENTS" : CODE})
164 #CODE=createGradientCode(f_interpolation, x, Q_r, loopindex=idx_full)
165 #insertCode(filename, { "SNIP_GRAD_REDUCED_ELEMENTS" : CODE})
166 #CODE=""
167 #for i in xrange(len(Q_faces)):
168 # CODE+="if (m_faceOffset[%s] > -1) {\n"%i
169 # CODE+=createGradientCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="m_faceOffset[%s]"%i)
170 # CODE+="\n} /* end of face %s */\n"%i
171 #insertCode(filename, { "SNIP_GRAD_FACES" : CODE})
172 #CODE=""
173 #for i in xrange(len(Q_r_faces)):
174 # CODE+="if (m_faceOffset[%s] > -1) {\n"%i
175 # CODE+=createGradientCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="m_faceOffset[%s]"%i)
176 # CODE+="\n} /* end of face %s */\n"%i
177 #insertCode(filename, { "SNIP_GRAD_REDUCED_FACES" : CODE})
178
179 #generate PDE assemblage
180 #CODE,PRECODE = makePDE(S, x, Q, W, DIM=DIM, system=True)
181 #insertCode(filename, { "SNIP_PDE_SYSTEM" : CODE , "SNIP_PDE_SYSTEM_PRE" : PRECODE } )
182 #CODE,PRECODE = makePDE(S, x, Q_r, W_r, DIM=DIM, system=True)
183 #insertCode(filename, { "SNIP_PDE_SYSTEM_REDUCED" : CODE , "SNIP_PDE_SYSTEM_REDUCED_PRE" : PRECODE })
184 #CODE,PRECODE = makePDE(S, x, Q, W, DIM=DIM, system=False)
185 #insertCode(filename, { "SNIP_PDE_SINGLE" : CODE , "SNIP_PDE_SINGLE_PRE" : PRECODE })
186 #CODE,PRECODE = makePDE(S, x, Q_r, W_r, DIM=DIM, system=False)
187 #insertCode(filename, { "SNIP_PDE_SINGLE_REDUCED" : CODE , "SNIP_PDE_SINGLE_REDUCED_PRE" : PRECODE })
188
189 #generate PDEBoundary assemblage
190 CODE,PRECODE = makePDEBC(S, x, Q_faces, W_faces, DIM=DIM, system=True)
191 insertCode(filename, extendDictionary( { "SNIP_PDEBC_SYSTEM_PRE" : PRECODE }, "SNIP_PDEBC_SYSTEM", CODE))
192 CODE,PRECODE = makePDEBC(S, x, Q_r_faces, W_r_faces, DIM=DIM, system=True)
193 insertCode(filename, extendDictionary( { "SNIP_PDEBC_SYSTEM_REDUCED_PRE" : PRECODE }, "SNIP_PDEBC_SYSTEM_REDUCED", CODE))
194 CODE,PRECODE = makePDEBC(S, x, Q_faces, W_faces, DIM=DIM, system=False)
195 insertCode(filename, extendDictionary( { "SNIP_PDEBC_SINGLE_PRE" : PRECODE }, "SNIP_PDEBC_SINGLE", CODE ))
196 CODE,PRECODE = makePDEBC(S, x, Q_r_faces, W_r_faces, DIM=DIM, system=False)
197 insertCode(filename, extendDictionary( { "SNIP_PDEBC_SINGLE_REDUCED_PRE" : PRECODE }, "SNIP_PDEBC_SINGLE_REDUCED", CODE))
198
199 def makePDEBC(S, x, Q, W, DIM, system=True):
200 GLOBAL_TMP={}
201 GLOBAL_N=0
202 PRECODE=""
203 CODE_OUT=[]
204
205 for n in range(len(Q)):
206 CODE="""
207 ///////////////
208 // process d //
209 ///////////////
210 if (add_EM_S) {
211 const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
212 """
213 if len(Q[n]) > 1:
214 CODE+="if (d.actsExpanded()) {\n"
215 if system: CODE+= """for (index_t k=0; k<numEq; k++) {
216 for (index_t m=0; m<numComp; m++) {
217 """
218 DATA_D=[]
219 CODE2=""
220 EM = {}
221 for i in range(len(S)):
222 for j in range(len(S)):
223 EM[(i,j)] = 0
224
225 for q in xrange(len(Q[n])):
226 A_name="d_%d"%(q,)
227 A=Symbol(A_name)
228 DATA_D.append(A)
229 if system:
230 CODE2+="const double %s = d_p[INDEX3(k, m, %s, numEq, numComp)];\n"%(A_name, q)
231 else:
232 CODE2+="const double %s = d_p[%s];\n"%(A_name, q)
233 for i in range(len(S)):
234 for j in range(len(S)):
235 EM[(i,j)] = EM[(i,j)] + (A * W[n][q] * S[j] * S[i]).subs( [ (x[jj], Q[n][q][jj]) for jj in xrange(DIM) ] )
236 CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP, system)
237
238 if system: CODE+="}\n }\n"
239 CODE+="} else { /* constant data */\n"
240 if system:
241 CODE+= """for (index_t k=0; k<numEq; k++) {
242 for (index_t m=0; m<numComp; m++) {
243 """
244 DATA_D=[]
245 CODE2=""
246 EM = {}
247 for i in range(len(S)):
248 for j in range(len(S)):
249 EM[(i,j)] = 0
250 if 1:
251 A_name="d_0"
252 A=Symbol(A_name)
253 DATA_D.append(A)
254 if system:
255 CODE2+="const double %s = d_p[INDEX2(k, m, numEq)];\n"%(A_name,)
256 else:
257 CODE2+="const double %s = d_p[0];\n"%(A_name,)
258 for q in xrange(len(Q[n])):
259 for i in range(len(S)):
260 for j in range(len(S)):
261 EM[(i,j)] = EM[(i,j)] + (A * W[n][q] * S[j] * S[i] ).subs( [ (x[jj], Q[n][q][jj]) for jj in xrange(DIM) ] )
262 CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP, system)
263 if system: CODE+="}\n}\n"
264 if len(Q[n]) > 1: CODE+="}\n"
265 CODE+="}\n"
266
267 #YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYy
268 CODE+="""
269 ///////////////
270 // process y //
271 ///////////////
272 if (add_EM_F) {
273 const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
274 """
275 if len(Q[n]) > 1:
276 CODE+="if (y.actsExpanded()) {\n"
277 if system: CODE+= "for (index_t k=0; k<numEq; k++) {\n"
278 DATA_Y=[]
279 CODE2=""
280 EM = {}
281 for i in range(len(S)):
282 EM[i] = 0
283
284 for q in xrange(len(Q[n])):
285 A_name="y_%d"%(q,)
286 A=Symbol(A_name)
287 DATA_Y.append(A)
288 if system:
289 CODE2+="const double %s = y_p[INDEX2(k, %s, numEq)];\n"%(A_name, q)
290 else:
291 CODE2+="const double %s = y_p[%s];\n"%(A_name, q)
292 for i in range(len(S)):
293 EM[i] = EM[i] + (A * W[n][q] * S[i]).subs( [ (x[jj], Q[n][q][jj]) for jj in xrange(DIM) ] )
294 CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP, system)
295
296 if system: CODE+="}\n"
297 CODE+="} else { /* constant data */\n"
298 if system:
299 CODE+= "for (index_t k=0; k<numEq; k++) {\n"
300 DATA_Y=[]
301 CODE2=""
302 EM = {}
303 for i in range(len(S)):
304 EM[i] = 0
305 if 1:
306 A_name="y_0"
307 A=Symbol(A_name)
308 DATA_Y.append(A)
309 if system:
310 CODE2+="const double %s = y_p[k];\n"%(A_name,)
311 else:
312 CODE2+="const double %s = y_p[0];\n"%(A_name,)
313 for q in xrange(len(Q[n])):
314 for i in range(len(S)):
315 EM[i] = EM[i] + (A * W[n][q] * S[i]).subs( [ (x[jj], Q[n][q][jj]) for jj in xrange(DIM) ] )
316 CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP, system)
317 if system: CODE+="}\n"
318 if len(Q[n]) > 1: CODE+="}\n"
319 CODE+="}\n"
320
321 CODE_OUT.append(CODE)
322 import operator
323 for k,v in sorted(GLOBAL_TMP.iteritems(), key=operator.itemgetter(1)):
324 PRECODE+="const double %s = %s;\n"%(v,ccode(k.evalf(n=DIGITS)))
325
326 return CODE_OUT, PRECODE
327
328 def extendDictionary(this, tag, items):
329 for i in range(len(items)):
330 this["%s_%s"%(tag,i)]=items[i]
331 return this
332
333 def insertCode(fn, replacement):
334 TOP_LINE_FMT="/* GENERATOR %s TOP"
335 BOTTOM_LINE_FMT="/* GENERATOR %s BOTTOM"
336 TABSIZE=4
337
338 # read code and create back up
339 bkp_stream=open(fn+".bkp",'w')
340 old_code=""
341 for l in open(fn, 'r').readlines():
342 s=l
343 if len(s)>0:
344 bkp_stream.write(s)
345 old_code+=s
346 bkp_stream.close()
347
348 for k in replacement.keys():
349 TOP_LINE=TOP_LINE_FMT%k
350 BOTTOM_LINE=BOTTOM_LINE_FMT%k
351 new_code=""
352 ignore=False
353 for l in old_code.split("\n"):
354 if l.strip().startswith(TOP_LINE):
355 ignore=True
356 N=len(l[:l.find(TOP_LINE)])
357 new_code+=l+"\n"
358 for rep_l in replacement[k].split("\n"):
359 #insert tabs:
360 if rep_l.find("/*")>-1 and rep_l.find("*/")>-1:
361 s=(rep_l[:rep_l.find("/*")] + rep_l[rep_l.find("*/")+2:]).strip()
362 else:
363 s=rep_l.strip()
364 if len(rep_l)>0:
365 if s.startswith("}"):
366 N=max(N-TABSIZE, 0)
367 if s.startswith("#"):
368 new_code+=rep_l + "\n"
369 else:
370 new_code+=" "*N + rep_l + "\n"
371 if s.endswith("{"):
372 N+=TABSIZE
373 elif l.strip().startswith(BOTTOM_LINE):
374 ignore=False
375 new_code+=l+"\n"
376 elif not ignore:
377 new_code+=l+"\n"
378 old_code=new_code
379 new_code=old_code[:len(old_code)-1]
380 open(fn,'w').write(new_code)
381 print "file %s updated."%fn
382
383 def optimizeEvaluate(F, x, Q):
384 F2=[]
385 for f in F:
386 f0=[]
387 for q in xrange(len(Q)):
388 tmp=f
389 for i in range(len(x)): tmp=tmp.subs(x[i], Q[q][i])
390 f0.append(tmp)
391 F2.append(f0)
392
393 # collect arguments:
394 args=[]
395 for f in F2:
396 for fq in f:
397 for s in fq.atoms():
398 if isinstance(s, Symbol):
399 if collect(fq, s, evaluate=False).has_key(s):
400 if s.name.startswith("f_"):
401 if args.count(s) == 0 and not (collect(fq, s, evaluate=False)[s]).is_zero: args.append(s)
402 c0={}
403 cc=0
404 F_new=[]
405 for f in F2:
406 fq_new=[]
407 for fq in f:
408 s1={}
409 for a in args:
410 s=collect(fq, a, evaluate=False)
411 if s.has_key(a):
412 k0=None
413 for k, v in c0.items():
414 if s[a] == v: k0=k
415 if k0==None:
416 k0=Symbol("tmp0_%s"%cc)
417 c0[k0]=s[a]
418 cc+=1
419 if not s1.has_key(k0): s1[k0]=0
420 s1[k0] = s1[k0] + a
421 tmp=0
422 for k,v in s1.items():
423 tmp = tmp + k * v
424 fq_new.append(tmp)
425 F_new.append(fq_new)
426
427 return F_new, args, c0
428
429 def createGradientCode(F, x, Q, gridoffset="", loopindex=(0,1,2)):
430 DIM=len(loopindex)
431 dF=[ expand(diff(F, x[i])) for i in range(DIM)]
432 dF, args, consts = optimizeEvaluate(dF, x, Q)
433 k=[]
434 N=[]
435 M1=""
436 TXT=""
437
438 for i in xrange(len(Q[0])):
439 if loopindex[i]>-1:
440 k.append("k%s"%i)
441 N.append("m_NE%s"%i)
442 if len(M1)==0:
443 M1+="m_N%s"%i
444 else:
445 if i < DIM-1: M1+=",m_N%s"%i
446
447 if len(k)>=2:
448 IDX2="INDEX%s(%s,%s)"%(len(k),",".join(k),",".join(N[:-1]))
449 else:
450 IDX2=k[0]
451
452 #for i in xrange(DIM-1,-1,-1):
453 # if loopindex[i] <= -1: TXT+="const index_t k%i = 0;\n"%i
454 for a,v in consts.items():
455 TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))
456 TXT+="#pragma omp parallel for"
457
458 for i in xrange(DIM-1,-1,-1):
459 if loopindex[i] > -1: TXT+="\nfor (index_t k%i=0; k%i < m_NE%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i])
460 for s in args:
461 k1=""
462 for i in xrange(DIM):
463 if s.name[2+i]=="0":
464 if loopindex[i]>-1:
465 k1+=",k%s"%i
466 elif loopindex[i] == -1:
467 k1+=",0"
468 elif loopindex[i] == -2:
469 k1+=",m_N%s-2"%i
470 else:
471 if loopindex[i]>-1:
472 k1+=",k%s+1"%i
473 elif loopindex[i] == -1:
474 k1+=",1"
475 elif loopindex[i] == -2:
476 k1+=",m_N%s-1"%i
477 TXT+="\nconst double* %s = in.getSampleDataRO(INDEX%s(%s, %s));"%(s.name,DIM,k1[1:],M1)
478 if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2
479 TXT+="\ndouble* o = out.getSampleDataRW(%s);"%IDX2
480 TXT+="\nfor (index_t i=0; i < numComp; ++i) {\n"
481 for q in xrange(len(Q)):
482 for i in range(DIM):
483 dFidx=dF[i][q]
484 for s in args:
485 dFidx=dFidx.subs(s.name, Symbol(s.name+"[i]"))
486 TXT+="o[INDEX3(i,%s,%s,numComp,%s)] = %s;\n"%(i,q,DIM,ccode(dFidx))
487 #TXT+="out[INDEX4(i,%s,%s,%s,numComp,%s,%s)] = %s;\n"%(i,q,IDX2,DIM,len(Q),ccode(dF[i][q]))
488 TXT+="} /* end of component loop i */\n"
489 for i in xrange(DIM):
490 if loopindex[i]>-1 : TXT+="} /* end of k%i loop */\n"%loopindex[i]
491 return TXT
492
493 def createCode(F, x, Q, gridoffset="", loopindex=(0,1,2)):
494 DIM=len(loopindex)
495 F, args, consts = optimizeEvaluate([F, ], x, Q)
496 F=F[0]
497 k=[]
498 N=[]
499 M1=""
500 TXT=""
501 for i in xrange(len(Q[0])):
502 if loopindex[i]>-1:
503 k.append("k%s"%i)
504 N.append("m_NE%s"%i)
505 if len(M1)==0:
506 M1+="m_N%s"%i
507 else:
508 if i < DIM-1: M1+=",m_N%s"%i
509
510 if len(k)>=2:
511 IDX2="INDEX%s(%s,%s)"%(len(k),",".join(k),",".join(N[:-1]))
512 else:
513 IDX2=k[0]
514
515 #for i in xrange(DIM-1,-1,-1):
516 # if loopindex[i] <= -1: TXT+="const index_t k%i = 0;\n"%i
517 for a,v in consts.items():
518 TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))
519 TXT+="#pragma omp parallel for"
520 for i in xrange(DIM-1,-1,-1):
521 if loopindex[i] > -1: TXT+="\nfor (index_t k%i=0; k%i < m_NE%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i])
522 for s in args:
523 k1=""
524 for i in xrange(DIM):
525 if s.name[2+i]=="0":
526 if loopindex[i]>-1:
527 k1+=",k%s"%i
528 elif loopindex[i] == -1:
529 k1+=",0"
530 elif loopindex[i] == -2:
531 k1+=",m_N%s-2"%i
532 else:
533 if loopindex[i]>-1:
534 k1+=",k%s+1"%i
535 elif loopindex[i] == -1:
536 k1+=",1"
537 elif loopindex[i] == -2:
538 k1+=",m_N%s-1"%i
539 TXT+="\nconst double* %s = in.getSampleDataRO(INDEX%s(%s, %s));"%(s.name,DIM,k1[1:],M1)
540 if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2
541 # for i in xrange(DIM):
542 # if loopindex[i]>-1: IDX2=gridoffset+"+k%s"%i
543 TXT+="\ndouble* o = out.getSampleDataRW(%s);"%IDX2
544 TXT+="\nfor (index_t i=0; i < numComp; ++i) {\n"
545 # interpolation to quadrature points
546 for q in xrange(len(Q)):
547 Fidx=F[q]
548 for s in args:
549 Fidx=Fidx.subs(s.name, Symbol(s.name+"[i]"))
550 TXT+="o[INDEX2(i,numComp,%s)] = %s;\n"%(q,ccode(Fidx))
551 TXT+="} /* end of component loop i */\n"
552 for i in xrange(DIM):
553 if loopindex[i]>-1 : TXT+="} /* end of k%i loop */\n"%loopindex[i]
554 return TXT
555
556 def createIntegrationCode(Q, W, gridoffset="", loopindex=(0,1,2)):
557 DIM=len(loopindex)
558 GLOBAL_TMP={}
559 LOCAL_TMP_E={}
560 TXT_E=""
561 for q in range(len(W)):
562 if not GLOBAL_TMP.has_key(W[q]):
563 GLOBAL_TMP[W[q]]=Symbol("w_%s"%len(GLOBAL_TMP))
564 A=Symbol("f_%s"%q)
565 TXT_E+="const double %s = in[INDEX3(i,%s,e, NCOMP,%s)];\n"%(A,q,len(W))
566 if not LOCAL_TMP_E.has_key(GLOBAL_TMP[W[q]]):
567 LOCAL_TMP_E[GLOBAL_TMP[W[q]]]=0
568 LOCAL_TMP_E[GLOBAL_TMP[W[q]]]=LOCAL_TMP_E[GLOBAL_TMP[W[q]]]+A
569 for p in LOCAL_TMP_E.keys():
570 TXT_E+="const double %s = %s;\n"%( p, ccode(LOCAL_TMP_E[p]))
571 print TXT_E
572 print GLOBAL_TMP
573 return ""
574 1/0
575 #=================
576
577
578
579 k3=""
580 for i in xrange(DIM-1,-1,-1):
581 if loopindex[i] > -1: k3+=",k%i"%loopindex[i]
582 TXT=""
583 for a,v in consts.items():
584 TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))
585 TXT+="#pragma omp parallel for private(i%s)"%k3
586 for i in xrange(DIM-1,-1,-1):
587 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])
588 TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"
589 print TXT
590 1/0
591 # interpolation to quadrature points
592 for q in xrange(len(Q)):
593 IDX2="INDEX%s(%s,%s)"%(DIM,k,N)
594 if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2
595 TXT+="out[INDEX3(i,%s,%s,NCOMP,%s)] = %s;\n"%(q,IDX2,len(Q),ccode(F[q]))
596 TXT+="} /* close component loop i */\n"
597 for i in xrange(DIM):
598 if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]
599 return TXT
600
601 def subPoint(g,x, p):
602 out=g
603 for i in range(len(x)): out=out.subs(x[i], p[i])
604 # return out.evalf()
605 return out
606
607 def generatePDECode(DATA_A, EM, GLOBAL_TMP, system=False):
608 LOCAL_TMP={}
609 LOCAL2_TMP={}
610 OUT=""
611 for p in EM.keys():
612 # first we collect the constant coefficients in the expression (outside the parallel region):
613 tmp={}
614 for a in DATA_A:
615 c= collect(EM[p], a, evaluate=False)
616 if c.has_key(a):
617 if not GLOBAL_TMP.has_key(c[a]):
618 GLOBAL_TMP[c[a]]=Symbol("w%s"%(len(GLOBAL_TMP),))
619 if tmp.has_key(GLOBAL_TMP[c[a]]):
620 tmp[GLOBAL_TMP[c[a]]]+=a
621 else:
622 tmp[GLOBAL_TMP[c[a]]]=a
623 # now we find temporary values which are combinations of the input data:
624 em_new=[]
625 for t in tmp.keys():
626 if tmp[t] in DATA_A:
627 tt = t * tmp[t]
628 else:
629 if not LOCAL_TMP.has_key(tmp[t]):
630 LOCAL_TMP[tmp[t]]=Symbol("tmp%s_0"%len(LOCAL_TMP))
631 tt= t * LOCAL_TMP[tmp[t]]
632 if not LOCAL2_TMP.has_key(tt):
633 LOCAL2_TMP[tt]=Symbol("tmp%s_1"%len(LOCAL2_TMP))
634 em_new.append(LOCAL2_TMP[tt])
635 EM[p]=em_new
636
637 for p in EM.keys():
638 sum_em=0
639 for t in EM[p]: sum_em=sum_em + t
640 if not isinstance(sum_em, int) :
641 if isinstance(p, int) and ( not isinstance(sum_em, int) ) :
642 if system:
643 OUT+="EM_F[INDEX2(k,%s,numEq)]+=%s;\n"%(p,ccode(sum_em))
644 else:
645 OUT+="EM_F[%s]+=%s;\n"%(p,ccode(sum_em))
646 else:
647 if system:
648 OUT+="EM_S[INDEX4(k,m,%s,%s,numEq,numComp,%s)]+=%s;\n"%(p[0],p[1],max([ qq[0] for qq in EM.keys()])+1,ccode(sum_em))
649 else:
650 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))
651
652 OUT2=""
653 for p in LOCAL_TMP:
654 OUT2+="const double %s = %s;\n"%(LOCAL_TMP[p],ccode(p))
655 for p in LOCAL2_TMP:
656 OUT2+="const double %s = %s;\n"%(LOCAL2_TMP[p],ccode(p))
657 return OUT2+OUT
658
659 def makePDE(S, x, Q, W, DIM=2, system=False):
660 GLOBAL_TMP={}
661 GLOBAL_N=0
662 PRECODE=""
663 CODE="""
664 ///////////////
665 // process A //
666 ///////////////
667 if (!A.isEmpty()) {
668 add_EM_S=true;
669 const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
670 """
671 if len(Q) > 1:
672 CODE+="if (A.actsExpanded()) {\n"
673 if system: CODE+= """for (index_t k=0; k<numEq; k++) {
674 for (index_t m=0; m<numComp; m++) {
675 """
676 DATA_A=[]
677 CODE2=""
678 EM = {}
679 for i in range(len(S)):
680 for j in range(len(S)):
681 EM[(i,j)] = 0
682
683 for q in xrange(len(Q)):
684 for di in range(DIM):
685 for dj in range(DIM):
686 A_name="A_%d%d_%d"%(di,dj,q)
687 A=Symbol(A_name)
688 DATA_A.append(A)
689 if system:
690 CODE2+="const double %s = A_p[INDEX5(k,%s,m,%s,%s, numEq,%s,numComp,%s)];\n"%(A_name, di, dj, q, DIM, DIM)
691 else:
692 CODE2+="const double %s = A_p[INDEX3(%s,%s,%s,%s,%s)];\n"%(A_name, di, dj, q, DIM, DIM)
693 for i in range(len(S)):
694 for j in range(len(S)):
695 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) ] )
696 CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP, system)
697
698 if system: CODE+="}\n}\n"
699 CODE+="} else { /* constant data */\n"
700 if system:
701 CODE+= """for (index_t k=0; k<numEq; k++) {
702 for (index_t m=0; m<numComp; m++) {
703 """
704 DATA_A=[]
705 CODE2=""
706 EM = {}
707 for i in range(len(S)):
708 for j in range(len(S)):
709 EM[(i,j)] = 0
710 for di in range(DIM):
711 for dj in range(DIM):
712 A_name="A_%d%d"%(di,dj)
713 A=Symbol(A_name)
714 DATA_A.append(A)
715 if system:
716 CODE2+="const double %s = A_p[INDEX4(k,%s,m,%s, numEq,%s, numComp)];\n"%(A_name, di, dj, DIM)
717 else:
718 CODE2+="const double %s = A_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, dj, DIM)
719 for q in xrange(len(Q)):
720 for i in range(len(S)):
721 for j in range(len(S)):
722 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) ] )
723 CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP, system)
724 if system: CODE+="}\n}\n"
725 if len(Q) > 1: CODE+="}\n"
726 CODE+="}\n"
727 #BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
728 CODE+="""
729 ///////////////
730 // process B //
731 ///////////////
732 if (!B.isEmpty()) {
733 add_EM_S=true;
734 const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
735 """
736 if len(Q) > 1:
737 CODE+="if (B.actsExpanded()) {\n"
738 if system: CODE+= """for (index_t k=0; k<numEq; k++) {
739 for (index_t m=0; m<numComp; m++) {
740 """
741 DATA_B=[]
742 CODE2=""
743 EM = {}
744 for i in range(len(S)):
745 for j in range(len(S)):
746 EM[(i,j)] = 0
747
748 for q in xrange(len(Q)):
749 for di in range(DIM):
750 A_name="B_%d_%d"%(di,q)
751 A=Symbol(A_name)
752 DATA_B.append(A)
753 if system:
754 CODE2+="const double %s = B_p[INDEX4(k,%s,m,%s, numEq,%s,numComp)];\n"%(A_name, di, q, DIM)
755 else:
756 CODE2+="const double %s = B_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, q, DIM)
757 for i in range(len(S)):
758 for j in range(len(S)):
759 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) ] )
760 CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP, system)
761
762 if system: CODE+="}\n}\n"
763 CODE+="} else { /* constant data */\n"
764 if system:
765 CODE+= """for (index_t k=0; k<numEq; k++) {
766 for (index_t m=0; m<numComp; m++) {
767 """
768 DATA_B=[]
769 CODE2=""
770 EM = {}
771 for i in range(len(S)):
772 for j in range(len(S)):
773 EM[(i,j)] = 0
774 for di in range(DIM):
775 A_name="B_%d"%(di)
776 A=Symbol(A_name)
777 DATA_B.append(A)
778 if system:
779 CODE2+="const double %s = B_p[INDEX3(k,%s,m, numEq, %s)];\n"%(A_name, di, DIM)
780 else:
781 CODE2+="const double %s = B_p[%s];\n"%(A_name, di)
782 for q in xrange(len(Q)):
783 for i in range(len(S)):
784 for j in range(len(S)):
785 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) ] )
786 CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP, system)
787 if system: CODE+="}\n}\n"
788 if len(Q) > 1: CODE+="}\n"
789 CODE+="}\n"
790 #CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
791 CODE+="""
792 ///////////////
793 // process C //
794 ///////////////
795 if (!C.isEmpty()) {
796 add_EM_S=true;
797 const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
798 """
799 if len(Q) > 1:
800 CODE+="if (C.actsExpanded()) {\n"
801 if system: CODE+= """for (index_t k=0; k<numEq; k++) {
802 for (index_t m=0; m<numComp; m++) {
803 """
804 DATA_C=[]
805 CODE2=""
806 EM = {}
807 for i in range(len(S)):
808 for j in range(len(S)):
809 EM[(i,j)] = 0
810
811 for q in xrange(len(Q)):
812 for dj in range(DIM):
813 A_name="C_%d_%d"%(dj,q)
814 A=Symbol(A_name)
815 DATA_C.append(A)
816 if system:
817 CODE2+="const double %s = C_p[INDEX4(k,m,%s, %s, numEq,numComp,%s)];\n"%(A_name, dj, q, DIM)
818 else:
819 CODE2+="const double %s = C_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj, q, DIM)
820 for i in range(len(S)):
821 for j in range(len(S)):
822 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) ] )
823 CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP, system)
824
825 if system: CODE+="}\n}\n"
826 CODE+="} else { /* constant data */\n"
827 if system:
828 CODE+= """for (index_t k=0; k<numEq; k++) {
829 for (index_t m=0; m<numComp; m++) {
830 """
831 DATA_C=[]
832 CODE2=""
833 EM = {}
834 for i in range(len(S)):
835 for j in range(len(S)):
836 EM[(i,j)] = 0
837 for dj in range(DIM):
838 A_name="C_%d"%(dj)
839 A=Symbol(A_name)
840 DATA_C.append(A)
841 if system:
842 CODE2+="const double %s = C_p[INDEX3(k, m, %s, numEq, numComp)];\n"%(A_name, dj)
843 else:
844 CODE2+="const double %s = C_p[%s];\n"%(A_name, dj)
845 for q in xrange(len(Q)):
846 for i in range(len(S)):
847 for j in range(len(S)):
848 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) ] )
849 CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP, system)
850 if system: CODE+="}\n}\n"
851 if len(Q) > 1: CODE+="}\n"
852 CODE+="}\n"
853 #DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD
854 CODE+="""
855 ///////////////
856 // process D //
857 ///////////////
858 if (!D.isEmpty()) {
859 add_EM_S=true;
860 const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
861 """
862 if len(Q) > 1:
863 CODE+="if (D.actsExpanded()) {\n"
864 if system: CODE+= """for (index_t k=0; k<numEq; k++) {
865 for (index_t m=0; m<numComp; m++) {
866 """
867 DATA_D=[]
868 CODE2=""
869 EM = {}
870 for i in range(len(S)):
871 for j in range(len(S)):
872 EM[(i,j)] = 0
873
874 for q in xrange(len(Q)):
875 A_name="D_%d"%(q,)
876 A=Symbol(A_name)
877 DATA_D.append(A)
878 if system:
879 CODE2+="const double %s = D_p[INDEX3(k, m, %s, numEq, numComp)];\n"%(A_name, q)
880 else:
881 CODE2+="const double %s = D_p[%s];\n"%(A_name, q)
882 for i in range(len(S)):
883 for j in range(len(S)):
884 EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )
885 CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP, system)
886
887 if system: CODE+="}\n }\n"
888 CODE+="} else { /* constant data */\n"
889 if system:
890 CODE+= """for (index_t k=0; k<numEq; k++) {
891 for (index_t m=0; m<numComp; m++) {
892 """
893 DATA_D=[]
894 CODE2=""
895 EM = {}
896 for i in range(len(S)):
897 for j in range(len(S)):
898 EM[(i,j)] = 0
899 if 1:
900 A_name="D_0"
901 A=Symbol(A_name)
902 DATA_D.append(A)
903 if system:
904 CODE2+="const double %s = D_p[INDEX2(k, m, numEq)];\n"%(A_name,)
905 else:
906 CODE2+="const double %s = D_p[0];\n"%(A_name,)
907 for q in xrange(len(Q)):
908 for i in range(len(S)):
909 for j in range(len(S)):
910 EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i] ).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )
911 CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP, system)
912 if system: CODE+="}\n}\n"
913 if len(Q) > 1: CODE+="}\n"
914 CODE+="}\n"
915
916 #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
917 CODE+="""
918 ///////////////
919 // process X //
920 ///////////////
921 if (!X.isEmpty()) {
922 add_EM_F=true;
923 const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
924 """
925 if len(Q) > 1:
926 CODE+="if (X.actsExpanded()) {\n"
927 if system: CODE+= "for (index_t k=0; k<numEq; k++) {\n"
928 DATA_X=[]
929 CODE2=""
930 EM = {}
931 for i in range(len(S)):
932 EM[i] = 0
933
934 for q in xrange(len(Q)):
935 for dj in range(DIM):
936 A_name="X_%d_%d"%(dj,q)
937 A=Symbol(A_name)
938 DATA_X.append(A)
939 if system:
940 CODE2+="const double %s = X_p[INDEX3(k, %s, %s, numEq, %s)];\n"%(A_name, dj, q, DIM)
941 else:
942 CODE2+="const double %s = X_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj,q,DIM)
943 for j in range(len(S)):
944 EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )
945 CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP, system)
946
947 if system: CODE+="}\n"
948 CODE+="} else { /* constant data */\n"
949 if system:
950 CODE+= "for (index_t k=0; k<numEq; k++) {\n"
951 DATA_X=[]
952 CODE2=""
953 EM = {}
954 for i in range(len(S)):
955 EM[i] = 0
956 for dj in range(DIM):
957 A_name="X_%d"%(dj)
958 A=Symbol(A_name)
959 DATA_X.append(A)
960 if system:
961 CODE2+="const double %s = X_p[INDEX2(k, %s, numEq)];\n"%(A_name, dj)
962 else:
963 CODE2+="const double %s = X_p[%s];\n"%(A_name, dj)
964 for q in xrange(len(Q)):
965 for j in range(len(S)):
966 EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )
967 CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP, system)
968 if system: CODE+="}\n"
969 if len(Q) > 1: CODE+="}\n"
970 CODE+="}\n"
971
972 #YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYy
973 CODE+="""
974 ///////////////
975 // process Y //
976 ///////////////
977 if (!Y.isEmpty()) {
978 add_EM_F=true;
979 const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
980 """
981 if len(Q) > 1:
982 CODE+="if (Y.actsExpanded()) {\n"
983 if system: CODE+= "for (index_t k=0; k<numEq; k++) {\n"
984 DATA_Y=[]
985 CODE2=""
986 EM = {}
987 for i in range(len(S)):
988 EM[i] = 0
989
990 for q in xrange(len(Q)):
991 A_name="Y_%d"%(q,)
992 A=Symbol(A_name)
993 DATA_Y.append(A)
994 if system:
995 CODE2+="const double %s = Y_p[INDEX2(k, %s, numEq)];\n"%(A_name, q)
996 else:
997 CODE2+="const double %s = Y_p[%s];\n"%(A_name, q)
998 for i in range(len(S)):
999 EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )
1000 CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP, system)
1001
1002 if system: CODE+="}\n"
1003 CODE+="} else { /* constant data */\n"
1004 if system:
1005 CODE+= "for (index_t k=0; k<numEq; k++) {\n"
1006 DATA_Y=[]
1007 CODE2=""
1008 EM = {}
1009 for i in range(len(S)):
1010 EM[i] = 0
1011 if 1:
1012 A_name="Y_0"
1013 A=Symbol(A_name)
1014 DATA_Y.append(A)
1015 if system:
1016 CODE2+="const double %s = Y_p[k];\n"%(A_name,)
1017 else:
1018 CODE2+="const double %s = Y_p[0];\n"%(A_name,)
1019 for q in xrange(len(Q)):
1020 for i in range(len(S)):
1021 EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )
1022 CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP, system)
1023 if system: CODE+="}\n"
1024 if len(Q) > 1: CODE+="}\n"
1025 CODE+="}\n"
1026
1027 import operator
1028 for k,v in sorted(GLOBAL_TMP.iteritems(), key=operator.itemgetter(1)):
1029 PRECODE+="const double %s = %s;\n"%(v,ccode(k.evalf(n=DIGITS)))
1030 return CODE, PRECODE
1031
1032 #filenames={2:"Rectangle.cpp", 3:"Brick.cpp"}
1033 filenames={3:"Brick.cpp"}
1034 for d in filenames.keys():
1035 generate(d, filenames[d])
1036

  ViewVC Help
Powered by ViewVC 1.1.26