/[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 4370 - (show annotations)
Fri Apr 19 06:15:24 2013 UTC (6 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 40668 byte(s)
WIP: ripley hand-optimisation & further rules to generator
-> drastically reduced number of constants
-> compile time of Brick.cpp less than half of what it was on savanna
-> additional runtime savings
-> to be continued...

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

  ViewVC Help
Powered by ViewVC 1.1.26