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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 6 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
2 ##############################################################################
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 from __future__ import print_function, division
18
19
20 """
21 this script generates the assemblage routine for the ripley rectangular grid solver
22
23
24 """
25 from sympy import *
26
27 DIGITS=20
28
29 #
30 # quadrature form in [0,1]
31 #
32
33 q0=1/sqrt(RealNumber(3))
34 qD1=[ (1-q0)/2, (1+q0)/2 ]
35 wD1=[RealNumber(1)/2, RealNumber(1)/2 ]
36
37
38 q0=1/sqrt(RealNumber(3))
39 qD1=[ (1-q0)/2, (1+q0)/2 ]
40 wD1=[RealNumber(1)/2, RealNumber(1)/2 ]
41
42 qD1_r=[ RealNumber(1)/2 ]
43 wD1_r= [RealNumber(1)]
44
45 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 h=[Symbol('h%s'%i) for i in range(DIM) ]
53 x=[Symbol('x%s'%i) for i in range(DIM) ]
54 # 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 for i in range(DIM):
64 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 for j in range(len(qD1)):
72 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 for j in range(len(qD1_r)):
78 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 for face in range(DIM):
89 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 for i in range(DIM):
95 # 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 for j in range(len(qD1)):
110 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 for i in range(DIM):
125 # 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 for j in range(len(qD1_r)):
135 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 for i in range(len(GridOffsetText)):
147 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 for i in range(len(Q_r_faces)):
156 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 for i in range(len(Q_r_faces)):
161 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 CODE+="\n} /* end of out_data_type branching */\n"
165 insertCode("Assemble_Interpolation_%sD.c"%DIM, { "SNIP" : CODE})
166
167 # gradient to Quad points
168 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 for i in range(len(Q_r_faces)):
174 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 for i in range(len(Q_r_faces)):
179 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 CODE+="\n} /* end of out_data_type branching */\n"
183 insertCode("Assemble_Gradient_%sD.c"%DIM, { "SNIP" : CODE})
184
185 #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
197
198 def insertCode(fn, replacement):
199 TOP_LINE_FMT="/* GENERATOR %s TOP"
200 BOTTOM_LINE_FMT="/* GENERATOR %s BOTTOM"
201 TAG=" "
202
203 # read code and create back up
204 bkp_stream=open(fn+".bkp",'w')
205 old_code=""
206 for l in open(fn, 'r').readlines():
207 s=l.strip()
208 if len(s)>0:
209 bkp_stream.write(s+"\n")
210 old_code+=s+"\n"
211 bkp_stream.close()
212
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 for l in old_code.split("\n"):
219 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 N=0
232 new_code=""
233 for l in old_code.split("\n"):
234 s=l.strip()
235 if s.find("/*")>-1 and s.find("*/")>-1 :
236 s2=(s[:s.find("/*")] + s[s.find("*/")+2:]).strip()
237 else:
238 s2=s
239 if len(s)>0:
240 if s2.startswith("}"):
241 N=max(N-1, 0)
242 new_code+=TAG*N + s + "\n"
243 if s2.endswith("{"):
244 N+=1
245 open(fn,'w').write(new_code)
246 print "file %s generated."%fn
247
248
249
250 def optimizeEvaluate(F, x, Q):
251 F2=[]
252 for f in F:
253 f0=[]
254 for q in range(len(Q)):
255 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 for i in range(len(Q[0])):
304 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 for i in range(DIM-1,-1,-1):
315 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 for i in range(DIM-1,-1,-1):
322 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 for i in range(DIM):
327 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 TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)
342 # interpolation to quadrature points
343 for q in range(len(Q)):
344 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 for i in range(DIM):
350 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 for i in range(len(Q[0])):
361 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 for i in range(DIM-1,-1,-1):
372 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 for i in range(DIM-1,-1,-1):
378 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 for i in range(DIM):
383 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 TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)
398 # interpolation to quadrature points
399 for q in range(len(Q)):
400 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 for i in range(DIM):
405 if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]
406 return TXT
407
408 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 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 # now we find temporary values which are combinations of the input data:
431 em_new=[]
432 for t in tmp.keys():
433 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 tt= t * LOCAL_TMP[tmp[t]]
439 if not LOCAL2_TMP.has_key(tt):
440 LOCAL2_TMP[tt]=Symbol("tmp%s_1"%len(LOCAL2_TMP))
441 em_new.append(LOCAL2_TMP[tt])
442 EM[p]=em_new
443
444 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 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 else:
453 if system:
454 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
458 OUT2=""
459 for p in LOCAL_TMP:
460 OUT2+=" const register double %s = %s;\n"%(LOCAL_TMP[p],ccode(p))
461 for p in LOCAL2_TMP:
462 OUT2+=" const register double %s = %s;\n"%(LOCAL2_TMP[p],ccode(p))
463 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 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 for dj in range(DIM):
491 A_name="A_%d%d_%d"%(di,dj,q)
492 A=Symbol(A_name)
493 DATA_A.append(A)
494 if system:
495 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 CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system)
502
503 if system: CODE+="}\n }\n"
504 CODE+="} else { /* constant data */\n"
505 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 EM[(i,j)] = 0
515 for di in range(DIM):
516 for dj in range(DIM):
517 A_name="A_%d%d"%(di,dj)
518 A=Symbol(A_name)
519 DATA_A.append(A)
520 if system:
521 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 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 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 A=Symbol(A_name)
556 DATA_B.append(A)
557 if system:
558 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 CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system)
565
566 if system: CODE+="}\n }\n"
567 CODE+="} else { /* constant data */\n"
568 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 EM[(i,j)] = 0
578 for di in range(DIM):
579 A_name="B_%d"%(di)
580 A=Symbol(A_name)
581 DATA_B.append(A)
582 if system:
583 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 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 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 A=Symbol(A_name)
618 DATA_C.append(A)
619 if system:
620 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 CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system)
627
628 if system: CODE+="}\n }\n"
629 CODE+="} else { /* constant data */\n"
630 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 EM[(i,j)] = 0
640 for dj in range(DIM):
641 A_name="C_%d"%(dj)
642 A=Symbol(A_name)
643 DATA_C.append(A)
644 if system:
645 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 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 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 A=Symbol(A_name)
679 DATA_D.append(A)
680 if system:
681 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 CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system)
688
689 if system: CODE+="}\n }\n"
690 CODE+="} else { /* constant data */\n"
691 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 EM[(i,j)] = 0
701 if 1:
702 A_name="D_0"
703 A=Symbol(A_name)
704 DATA_D.append(A)
705 if system:
706 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 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 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 A=Symbol(A_name)
740 DATA_X.append(A)
741 if system:
742 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 CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system)
748
749 if system: CODE+="}\n"
750 CODE+="} else { /* constant data */\n"
751 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 EM[i] = 0
759 for dj in range(DIM):
760 A_name="X_%d"%(dj)
761 A=Symbol(A_name)
762 DATA_X.append(A)
763 if system:
764 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 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 EM[i] = 0
792
793 for q in range(len(Q)):
794 A_name="Y_%d"%(q,)
795 A=Symbol(A_name)
796 DATA_Y.append(A)
797 if system:
798 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 CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system)
804
805 if system: CODE+="}\n"
806 CODE+="} else { /* constant data */\n"
807 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 EM[i] = 0
815 if 1:
816 A_name="Y_0"
817 A=Symbol(A_name)
818 DATA_Y.append(A)
819 if system:
820 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 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 for d in [3]:
836 generate(d)

  ViewVC Help
Powered by ViewVC 1.1.26