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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4575 by caltinay, Wed Feb 1 06:16:25 2012 UTC revision 4576 by sshaw, Mon Dec 9 23:35:30 2013 UTC
# Line 30  print "1D reduced quadrature weights =", Line 30  print "1D reduced quadrature weights =",
30    
31  def generate(DIM):  def generate(DIM):
32     # variables     # variables
33     h=[Symbol('h%s'%i) for i in xrange(DIM) ]     h=[Symbol('h%s'%i) for i in range(DIM) ]
34     x=[Symbol('x%s'%i) for i in xrange(DIM) ]     x=[Symbol('x%s'%i) for i in range(DIM) ]
35     #  shape functions: labeling differently from finley!!!!     #  shape functions: labeling differently from finley!!!!
36     S=[ 1 ]     S=[ 1 ]
37     GridOffset=[[]]     GridOffset=[[]]
# Line 41  def generate(DIM): Line 41  def generate(DIM):
41     W = [ 1 ]     W = [ 1 ]
42     Q_r = [ [] ]     Q_r = [ [] ]
43     W_r = [ 1 ]     W_r = [ 1 ]
44     for i in xrange(DIM):     for i in range(DIM):
45         idx_full.append(i)         idx_full.append(i)
46         S = [ s * (1-x[i]/h[i])  for s in S] + [ s * x[i]/h[i]  for s in S]         S = [ s * (1-x[i]/h[i])  for s in S] + [ s * x[i]/h[i]  for s in S]
47         GridOffset = [ s + [0]  for s in GridOffset] + [ s + [1]  for s in GridOffset]         GridOffset = [ s + [0]  for s in GridOffset] + [ s + [1]  for s in GridOffset]
# Line 49  def generate(DIM): Line 49  def generate(DIM):
49    
50         # generate quadrature points in element         # generate quadrature points in element
51         Qnew, Wnew =[], []         Qnew, Wnew =[], []
52         for j in xrange(len(qD1)):         for j in range(len(qD1)):
53                Qnew += [ q + [qD1[j]*h[i],] for q in Q ]                Qnew += [ q + [qD1[j]*h[i],] for q in Q ]
54                Wnew += [ w * wD1[j]*h[i] for w in W ]                Wnew += [ w * wD1[j]*h[i] for w in W ]
55         Q, W=Qnew, Wnew         Q, W=Qnew, Wnew
56    
57         Qnew, Wnew =[], []         Qnew, Wnew =[], []
58         for j in xrange(len(qD1_r)):         for j in range(len(qD1_r)):
59                 Qnew += [ q + [qD1_r[j]*h[i],] for q in Q_r ]                 Qnew += [ q + [qD1_r[j]*h[i],] for q in Q_r ]
60                 Wnew += [ w * wD1_r[j]*h[i] for w in W_r ]                 Wnew += [ w * wD1_r[j]*h[i] for w in W_r ]
61         Q_r, W_r=Qnew, Wnew         Q_r, W_r=Qnew, Wnew
# Line 66  def generate(DIM): Line 66  def generate(DIM):
66     Q_r_faces=[]     Q_r_faces=[]
67     W_r_faces=[]     W_r_faces=[]
68     idx_faces=[]     idx_faces=[]
69     for face in xrange(DIM):     for face in range(DIM):
70         Q_left, W_left = [ [] ], [ 1 ]         Q_left, W_left = [ [] ], [ 1 ]
71         Q_left_r, W_left_r = [ [] ], [ 1 ]         Q_left_r, W_left_r = [ [] ], [ 1 ]
72         Q_right, W_right = [ [] ], [ 1 ]         Q_right, W_right = [ [] ], [ 1 ]
73         Q_right_r, W_right_r = [ [] ], [ 1 ]         Q_right_r, W_right_r = [ [] ], [ 1 ]
74         idx_left,idx_right=[],[]         idx_left,idx_right=[],[]
75         for i in xrange(DIM):           for i in range(DIM):  
76            # generate quadrature points in element            # generate quadrature points in element
77            Q_left_new, W_left_new =[], []            Q_left_new, W_left_new =[], []
78            Q_right_new, W_right_new =[], []            Q_right_new, W_right_new =[], []
# Line 87  def generate(DIM): Line 87  def generate(DIM):
87            else:            else:
88               idx_left.append(i)               idx_left.append(i)
89               idx_right.append(i)               idx_right.append(i)
90               for j in xrange(len(qD1)):               for j in range(len(qD1)):
91                   Q_left_new += [ q + [qD1[j]*h[i],] for q in Q_left ]                   Q_left_new += [ q + [qD1[j]*h[i],] for q in Q_left ]
92                   W_left_new += [ w * wD1[j]*h[i] for w in W_left ]                   W_left_new += [ w * wD1[j]*h[i] for w in W_left ]
93                   Q_right_new += [ q + [qD1[j]*h[i],] for q in Q_right ]                   Q_right_new += [ q + [qD1[j]*h[i],] for q in Q_right ]
# Line 102  def generate(DIM): Line 102  def generate(DIM):
102         Q_left_r, W_left_r = [ [] ], [ 1 ]         Q_left_r, W_left_r = [ [] ], [ 1 ]
103         Q_right, W_right = [ [] ], [ 1 ]         Q_right, W_right = [ [] ], [ 1 ]
104         Q_right_r, W_right_r = [ [] ], [ 1 ]         Q_right_r, W_right_r = [ [] ], [ 1 ]
105         for i in xrange(DIM):         for i in range(DIM):
106            # generate quadrature points in element            # generate quadrature points in element
107            Q_left_new, W_left_new =[], []            Q_left_new, W_left_new =[], []
108            Q_right_new, W_right_new =[], []            Q_right_new, W_right_new =[], []
# Line 112  def generate(DIM): Line 112  def generate(DIM):
112               Q_right_new += [ q + [ h[i], ] for q in Q_right ]               Q_right_new += [ q + [ h[i], ] for q in Q_right ]
113               W_right_new += [ w * 1. for w in W_right ]               W_right_new += [ w * 1. for w in W_right ]
114            else:            else:
115               for j in xrange(len(qD1_r)):               for j in range(len(qD1_r)):
116                   Q_left_new += [ q + [qD1_r[j]*h[i],] for q in Q_left ]                   Q_left_new += [ q + [qD1_r[j]*h[i],] for q in Q_left ]
117                   W_left_new += [ w * wD1_r[j]*h[i] for w in W_left ]                   W_left_new += [ w * wD1_r[j]*h[i] for w in W_left ]
118                   Q_right_new += [ q + [qD1_r[j]*h[i],] for q in Q_right ]                   Q_right_new += [ q + [qD1_r[j]*h[i],] for q in Q_right ]
# Line 124  def generate(DIM): Line 124  def generate(DIM):
124                
125     # the interpolation function     # the interpolation function
126     f_interpolation = 0     f_interpolation = 0
127     for i in xrange(len(GridOffsetText)):     for i in range(len(GridOffsetText)):
128       f_interpolation = f_interpolation + S[i] * Symbol("f_%s"%GridOffsetText[i])       f_interpolation = f_interpolation + S[i] * Symbol("f_%s"%GridOffsetText[i])
129    
130     # interpolation to Quad points     # interpolation to Quad points
# Line 133  def generate(DIM): Line 133  def generate(DIM):
133     CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"     CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"
134     CODE+=createCode(f_interpolation, x, Q_r, loopindex=idx_full)     CODE+=createCode(f_interpolation, x, Q_r, loopindex=idx_full)
135     CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"     CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"
136     for i in xrange(len(Q_r_faces)):     for i in range(len(Q_r_faces)):
137          CODE+="if (face_offset(%s)>-1) {\n"%i          CODE+="if (face_offset(%s)>-1) {\n"%i
138          CODE+=createCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)          CODE+=createCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)
139          CODE+="\n} /* end of face %s */\n"%i          CODE+="\n} /* end of face %s */\n"%i
140     CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"     CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"
141     for i in xrange(len(Q_r_faces)):     for i in range(len(Q_r_faces)):
142          CODE+="if (face_offset(%s)>-1) {\n"%i          CODE+="if (face_offset(%s)>-1) {\n"%i
143          CODE+=createCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)          CODE+=createCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)
144          CODE+="\n} /* end of face %s */\n"%i          CODE+="\n} /* end of face %s */\n"%i
# Line 151  def generate(DIM): Line 151  def generate(DIM):
151     CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"     CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"
152     CODE+=createGradientCode(f_interpolation, x, Q_r, loopindex=idx_full)     CODE+=createGradientCode(f_interpolation, x, Q_r, loopindex=idx_full)
153     CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"     CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"
154     for i in xrange(len(Q_r_faces)):     for i in range(len(Q_r_faces)):
155          CODE+="if (face_offset(%s)>-1) {\n"%i          CODE+="if (face_offset(%s)>-1) {\n"%i
156          CODE+=createGradientCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)          CODE+=createGradientCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)
157          CODE+="\n} /* end of face %s */\n"%i          CODE+="\n} /* end of face %s */\n"%i
158     CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"     CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"
159     for i in xrange(len(Q_r_faces)):     for i in range(len(Q_r_faces)):
160          CODE+="if (face_offset(%s)>-1) {\n"%i          CODE+="if (face_offset(%s)>-1) {\n"%i
161          CODE+=createGradientCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)          CODE+=createGradientCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)
162          CODE+="\n} /* end of face %s */\n"%i          CODE+="\n} /* end of face %s */\n"%i
# Line 196  def insertCode(fn, replacement): Line 196  def insertCode(fn, replacement):
196          BOTTOM_LINE=BOTTOM_LINE_FMT%k          BOTTOM_LINE=BOTTOM_LINE_FMT%k
197          new_code=""          new_code=""
198          ignore=False          ignore=False
199          for l in old_code.split("\n"):              for l in old_code.split("\n"):          
200             s=l.strip()             s=l.strip()
201             if len(s)>0:             if len(s)>0:
202                if s.startswith(TOP_LINE):                if s.startswith(TOP_LINE):
# Line 214  def insertCode(fn, replacement): Line 214  def insertCode(fn, replacement):
214      for l in old_code.split("\n"):      for l in old_code.split("\n"):
215         s=l.strip()         s=l.strip()
216         if s.find("/*")>-1 and s.find("*/")>-1 :           if s.find("/*")>-1 and s.find("*/")>-1 :  
217          s2=(s[:s.find("/*")] + s[s.find("*/")+2:]).strip()              s2=(s[:s.find("/*")] + s[s.find("*/")+2:]).strip()
218         else:         else:
219          s2=s              s2=s
220         if len(s)>0:         if len(s)>0:
221            if s2.startswith("}"):            if s2.startswith("}"):
222               N=max(N-1, 0)               N=max(N-1, 0)
# Line 232  def optimizeEvaluate(F, x, Q): Line 232  def optimizeEvaluate(F, x, Q):
232    F2=[]    F2=[]
233    for f in F:    for f in F:
234       f0=[]       f0=[]
235       for q in xrange(len(Q)):       for q in range(len(Q)):
236          tmp=f          tmp=f
237          for i in range(len(x)): tmp=tmp.subs(x[i], Q[q][i])          for i in range(len(x)): tmp=tmp.subs(x[i], Q[q][i])
238          f0.append(tmp)          f0.append(tmp)
# Line 281  def createGradientCode(F, x, Q, gridoffs Line 281  def createGradientCode(F, x, Q, gridoffs
281     k=""     k=""
282     N=""     N=""
283     M1=""     M1=""
284     for i in xrange(len(Q[0])):     for i in range(len(Q[0])):
285       if len(k)==0:       if len(k)==0:
286           k+="k%s"%i           k+="k%s"%i
287           N+="N%s"%i           N+="N%s"%i
# Line 292  def createGradientCode(F, x, Q, gridoffs Line 292  def createGradientCode(F, x, Q, gridoffs
292           if i < DIM-1: M1+=",M%s"%i           if i < DIM-1: M1+=",M%s"%i
293        
294     k3=""     k3=""
295     for i in xrange(DIM-1,-1,-1):     for i in range(DIM-1,-1,-1):
296            if loopindex[i] > -1: k3+=",k%i"%loopindex[i]            if loopindex[i] > -1: k3+=",k%i"%loopindex[i]
297     TXT=""     TXT=""
298     for a,v in  consts.items():     for a,v in  consts.items():
299        TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))        TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))
300     TXT+="#pragma omp parallel for private(i%s)"%k3     TXT+="#pragma omp parallel for private(i%s)"%k3
301    
302     for i in xrange(DIM-1,-1,-1):     for i in range(DIM-1,-1,-1):
303            if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i])            if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i])
304     TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"     TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"
305     for s in args:     for s in args:
306              k1=""              k1=""
307              for i in xrange(DIM):              for i in range(DIM):
308                if s.name[2+i]=="0":                if s.name[2+i]=="0":
309                    if loopindex[i]>-1:                    if loopindex[i]>-1:
310                       k1+=",k%s"%i                       k1+=",k%s"%i
# Line 321  def createGradientCode(F, x, Q, gridoffs Line 321  def createGradientCode(F, x, Q, gridoffs
321                       k1+=",M%s-1"%i                       k1+=",M%s-1"%i
322              TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)              TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)
323     #  interpolation to quadrature points     #  interpolation to quadrature points
324     for q in xrange(len(Q)):     for q in range(len(Q)):
325           IDX2="INDEX%s(%s,%s)"%(DIM,k,N)           IDX2="INDEX%s(%s,%s)"%(DIM,k,N)
326           if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2           if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2
327           for i in range(DIM):           for i in range(DIM):
328               TXT+="out[INDEX4(i,%s,%s,%s,NCOMP,%s,%s)] = %s;\n"%(i,q,IDX2,DIM,len(Q),ccode(dF[i][q]))               TXT+="out[INDEX4(i,%s,%s,%s,NCOMP,%s,%s)] = %s;\n"%(i,q,IDX2,DIM,len(Q),ccode(dF[i][q]))
329     TXT+="} /* close component loop i */\n"     TXT+="} /* close component loop i */\n"
330     for i in xrange(DIM):     for i in range(DIM):
331             if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]             if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]
332     return TXT     return TXT
333    
# Line 338  def createCode(F, x, Q, gridoffset="", l Line 338  def createCode(F, x, Q, gridoffset="", l
338     k=""     k=""
339     N=""     N=""
340     M1=""     M1=""
341     for i in xrange(len(Q[0])):     for i in range(len(Q[0])):
342       if len(k)==0:       if len(k)==0:
343           k+="k%s"%i           k+="k%s"%i
344           N+="N%s"%i           N+="N%s"%i
# Line 349  def createCode(F, x, Q, gridoffset="", l Line 349  def createCode(F, x, Q, gridoffset="", l
349           if i < DIM-1: M1+=",M%s"%i           if i < DIM-1: M1+=",M%s"%i
350        
351     k3=""     k3=""
352     for i in xrange(DIM-1,-1,-1):     for i in range(DIM-1,-1,-1):
353            if loopindex[i] > -1: k3+=",k%i"%loopindex[i]            if loopindex[i] > -1: k3+=",k%i"%loopindex[i]
354     TXT=""     TXT=""
355     for a,v in  consts.items():     for a,v in  consts.items():
356        TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))        TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))
357     TXT+="#pragma omp parallel for private(i%s)"%k3     TXT+="#pragma omp parallel for private(i%s)"%k3
358     for i in xrange(DIM-1,-1,-1):     for i in range(DIM-1,-1,-1):
359            if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i])            if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i])
360     TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"     TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"
361     for s in args:     for s in args:
362              k1=""              k1=""
363              for i in xrange(DIM):              for i in range(DIM):
364                if s.name[2+i]=="0":                if s.name[2+i]=="0":
365                    if loopindex[i]>-1:                    if loopindex[i]>-1:
366                       k1+=",k%s"%i                       k1+=",k%s"%i
# Line 377  def createCode(F, x, Q, gridoffset="", l Line 377  def createCode(F, x, Q, gridoffset="", l
377                       k1+=",M%s-1"%i                       k1+=",M%s-1"%i
378              TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)              TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1)
379     #  interpolation to quadrature points     #  interpolation to quadrature points
380     for q in xrange(len(Q)):     for q in range(len(Q)):
381        IDX2="INDEX%s(%s,%s)"%(DIM,k,N)        IDX2="INDEX%s(%s,%s)"%(DIM,k,N)
382        if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2        if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2
383        TXT+="out[INDEX3(i,%s,%s,NCOMP,%s)] = %s;\n"%(q,IDX2,len(Q),ccode(F[q]))        TXT+="out[INDEX3(i,%s,%s,NCOMP,%s)] = %s;\n"%(q,IDX2,len(Q),ccode(F[q]))
384     TXT+="} /* close component loop i */\n"     TXT+="} /* close component loop i */\n"
385     for i in xrange(DIM):     for i in range(DIM):
386             if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]             if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]
387     return TXT     return TXT
388        
# Line 392  def subPoint(g,x, p): Line 392  def subPoint(g,x, p):
392       # return out.evalf()       # return out.evalf()
393       return out       return out
394    
395  def generatePDECode(DATA_A, EM,GLOBAL_TMP, system=False):    def generatePDECode(DATA_A, EM,GLOBAL_TMP, system=False):        
396      LOCAL_TMP={}          LOCAL_TMP={}
397      LOCAL2_TMP={}            LOCAL2_TMP={}
398      OUT=""          OUT=""
399      for p in EM.keys():          for p in EM.keys():
400          # first we collect the constant coefficients in the expression (outside the parallel region):              # first we collect the constant coefficients in the expression (outside the parallel region):
401          tmp={}              tmp={}
402          for a in DATA_A:              for a in DATA_A:
403              c=  collect(EM[p], a, evaluate=False)                  c=  collect(EM[p], a, evaluate=False)
404              if c.has_key(a):                  if c.has_key(a):
405                 if not GLOBAL_TMP.has_key(c[a]):                     if not GLOBAL_TMP.has_key(c[a]):
406                         GLOBAL_TMP[c[a]]=Symbol("w%s"%(len(GLOBAL_TMP),))                             GLOBAL_TMP[c[a]]=Symbol("w%s"%(len(GLOBAL_TMP),))
407             if tmp.has_key(GLOBAL_TMP[c[a]]):                     if tmp.has_key(GLOBAL_TMP[c[a]]):
408                     tmp[GLOBAL_TMP[c[a]]]+=a                         tmp[GLOBAL_TMP[c[a]]]+=a
409             else:                     else:
410                  tmp[GLOBAL_TMP[c[a]]]=a                          tmp[GLOBAL_TMP[c[a]]]=a
411              # now we find temporary values which are combinations of the input data:              # now we find temporary values which are combinations of the input data:
412              em_new=[]              em_new=[]
413              for t in tmp.keys():              for t in tmp.keys():
414              if tmp[t] in DATA_A:                  if tmp[t] in DATA_A:
415             tt = t * tmp[t]                     tt = t * tmp[t]
416          else:                  else:
417                 if not LOCAL_TMP.has_key(tmp[t]):                     if not LOCAL_TMP.has_key(tmp[t]):
418                LOCAL_TMP[tmp[t]]=Symbol("tmp%s_0"%len(LOCAL_TMP))                        LOCAL_TMP[tmp[t]]=Symbol("tmp%s_0"%len(LOCAL_TMP))
419                     tt= t * LOCAL_TMP[tmp[t]]                     tt= t * LOCAL_TMP[tmp[t]]
420                  if not LOCAL2_TMP.has_key(tt):                  if not LOCAL2_TMP.has_key(tt):
421                LOCAL2_TMP[tt]=Symbol("tmp%s_1"%len(LOCAL2_TMP))                        LOCAL2_TMP[tt]=Symbol("tmp%s_1"%len(LOCAL2_TMP))
422                  em_new.append(LOCAL2_TMP[tt])                  em_new.append(LOCAL2_TMP[tt])
423              EM[p]=em_new              EM[p]=em_new
424                                                
425          for p in EM.keys():          for p in EM.keys():
426              sum_em=0              sum_em=0
427              for t in EM[p]: sum_em=sum_em + t              for t in EM[p]: sum_em=sum_em + t
428              if isinstance(p, int) :              if isinstance(p, int) :
429            if system:                if system:
430               OUT+="  EM_F[INDEX2(k,%s,p.numEqu)]+=%s;\n"%(p,ccode(sum_em))                   OUT+="  EM_F[INDEX2(k,%s,p.numEqu)]+=%s;\n"%(p,ccode(sum_em))
431            else:                else:
432           OUT+="  EM_F[%s]+=%s;\n"%(p,ccode(sum_em))                   OUT+="  EM_F[%s]+=%s;\n"%(p,ccode(sum_em))
433              else:              else:
434                if system:                if system:
435               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))                   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))
436            else:                else:
437           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))                   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))
438    
439      OUT2=""          OUT2=""
440      for p in LOCAL_TMP:          for p in LOCAL_TMP:
441           OUT2+="  const register double %s = %s;\n"%(LOCAL_TMP[p],ccode(p))               OUT2+="  const register double %s = %s;\n"%(LOCAL_TMP[p],ccode(p))
442          for p in LOCAL2_TMP:          for p in LOCAL2_TMP:
443           OUT2+="  const register double %s = %s;\n"%(LOCAL2_TMP[p],ccode(p))               OUT2+="  const register double %s = %s;\n"%(LOCAL2_TMP[p],ccode(p))
444          return OUT2+OUT          return OUT2+OUT
445    
446  def makePDE(S, x, Q, W, DIM=2, system=False):  def makePDE(S, x, Q, W, DIM=2, system=False):
# Line 463  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 463  def makePDE(S, x, Q, W, DIM=2, system=Fa
463          CODE2=""          CODE2=""
464          EM = {}          EM = {}
465          for i in range(len(S)):          for i in range(len(S)):
466          for j in range(len(S)):              for j in range(len(S)):
467              EM[(i,j)] = 0                  EM[(i,j)] = 0
468                                
469          for q in xrange(len(Q)):          for q in range(len(Q)):
470          for di in range(DIM):              for di in range(DIM):
471                 for dj in range(DIM):                 for dj in range(DIM):
472                A_name="A_%d%d_%d"%(di,dj,q)                      A_name="A_%d%d_%d"%(di,dj,q)
473                    A=Symbol(A_name)                    A=Symbol(A_name)
474                    DATA_A.append(A)                    DATA_A.append(A)
475                    if system:                    if system:
476                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)                        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)
477            else:                    else:
478                CODE2+="const register double %s = A_p[INDEX3(%s,%s,%s,%s,%s)];\n"%(A_name, di, dj, q, DIM, DIM)                        CODE2+="const register double %s = A_p[INDEX3(%s,%s,%s,%s,%s)];\n"%(A_name, di, dj, q, DIM, DIM)
479            for i in range(len(S)):                     for i in range(len(S)):
480                    for j in range(len(S)):                        for j in range(len(S)):
481                        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) ] )                            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) ] )
482          CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system)          CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system)
483        
484      if system: CODE+="}\n }\n"          if system: CODE+="}\n }\n"
485      CODE+="} else { /* constant data */\n"          CODE+="} else { /* constant data */\n"
486     if system:     if system:
487         CODE+= """for (k=0;k<p.numEqu;k++) {         CODE+= """for (k=0;k<p.numEqu;k++) {
488                              for (m=0;m<p.numComp;m++) {                              for (m=0;m<p.numComp;m++) {
# Line 492  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 492  def makePDE(S, x, Q, W, DIM=2, system=Fa
492     EM = {}     EM = {}
493     for i in range(len(S)):     for i in range(len(S)):
494         for j in range(len(S)):         for j in range(len(S)):
495        EM[(i,j)] = 0            EM[(i,j)] = 0
496     for di in range(DIM):     for di in range(DIM):
497         for dj in range(DIM):         for dj in range(DIM):
498           A_name="A_%d%d"%(di,dj)                 A_name="A_%d%d"%(di,dj)
499               A=Symbol(A_name)               A=Symbol(A_name)
500               DATA_A.append(A)               DATA_A.append(A)
501               if system:               if system:
502                CODE2+="const register double %s = A_p[INDEX4(k,%s,m,%s p.numEqu,%s, p.numComp)];\n"%(A_name, di, dj, DIM)                        CODE2+="const register double %s = A_p[INDEX4(k,%s,m,%s p.numEqu,%s, p.numComp)];\n"%(A_name, di, dj, DIM)
503           else:               else:
504                CODE2+="const register double %s = A_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, dj, DIM)                        CODE2+="const register double %s = A_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, dj, DIM)
505               for q in xrange(len(Q)):               for q in range(len(Q)):
506            for i in range(len(S)):                     for i in range(len(S)):
507                    for j in range(len(S)):                        for j in range(len(S)):
508                        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) ] )                            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) ] )
509     CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system)         CODE+=CODE2+generatePDECode(DATA_A, EM, GLOBAL_TMP,system)    
510     if system:  CODE+="}\n }\n"     if system:  CODE+="}\n }\n"
511     if len(Q) > 1: CODE+="}\n"     if len(Q) > 1: CODE+="}\n"
# Line 527  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 527  def makePDE(S, x, Q, W, DIM=2, system=Fa
527          CODE2=""          CODE2=""
528          EM = {}          EM = {}
529          for i in range(len(S)):          for i in range(len(S)):
530          for j in range(len(S)):              for j in range(len(S)):
531              EM[(i,j)] = 0                  EM[(i,j)] = 0
532                                
533          for q in xrange(len(Q)):          for q in range(len(Q)):
534          for di in range(DIM):              for di in range(DIM):
535                A_name="B_%d_%d"%(di,q)                      A_name="B_%d_%d"%(di,q)
536                    A=Symbol(A_name)                    A=Symbol(A_name)
537                    DATA_B.append(A)                    DATA_B.append(A)
538                    if system:                    if system:
539                CODE2+="const register double %s = B_p[INDEX4(k,%s,m,%s, p.numEqu,%s,p.numComp)];\n"%(A_name, di,  q, DIM)                        CODE2+="const register double %s = B_p[INDEX4(k,%s,m,%s, p.numEqu,%s,p.numComp)];\n"%(A_name, di,  q, DIM)
540            else:                    else:
541                CODE2+="const register double %s = B_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, q, DIM)                        CODE2+="const register double %s = B_p[INDEX2(%s,%s,%s)];\n"%(A_name, di, q, DIM)
542            for i in range(len(S)):                     for i in range(len(S)):
543                    for j in range(len(S)):                        for j in range(len(S)):
544                        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) ] )                            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) ] )
545          CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system)          CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system)
546                    
547      if system: CODE+="}\n }\n"          if system: CODE+="}\n }\n"
548      CODE+="} else { /* constant data */\n"          CODE+="} else { /* constant data */\n"
549     if system:     if system:
550         CODE+= """for (k=0;k<p.numEqu;k++) {         CODE+= """for (k=0;k<p.numEqu;k++) {
551                              for (m=0;m<p.numComp;m++) {                              for (m=0;m<p.numComp;m++) {
# Line 555  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 555  def makePDE(S, x, Q, W, DIM=2, system=Fa
555     EM = {}     EM = {}
556     for i in range(len(S)):     for i in range(len(S)):
557         for j in range(len(S)):         for j in range(len(S)):
558        EM[(i,j)] = 0            EM[(i,j)] = 0
559     for di in range(DIM):     for di in range(DIM):
560           A_name="B_%d"%(di)                 A_name="B_%d"%(di)
561               A=Symbol(A_name)               A=Symbol(A_name)
562               DATA_B.append(A)               DATA_B.append(A)
563               if system:               if system:
564                CODE2+="const register double %s = B_p[INDEX3(k,%s,m, p.numEqu,%s)];\n"%(A_name, di,  DIM)                        CODE2+="const register double %s = B_p[INDEX3(k,%s,m, p.numEqu,%s)];\n"%(A_name, di,  DIM)
565           else:               else:
566                CODE2+="const register double %s = B_p[%s];\n"%(A_name, di)                        CODE2+="const register double %s = B_p[%s];\n"%(A_name, di)
567               for q in xrange(len(Q)):               for q in range(len(Q)):
568            for i in range(len(S)):                     for i in range(len(S)):
569                    for j in range(len(S)):                        for j in range(len(S)):
570                        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) ] )                            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) ] )
571     CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system)         CODE+=CODE2+generatePDECode(DATA_B, EM, GLOBAL_TMP,system)    
572     if system:  CODE+="}\n }\n"     if system:  CODE+="}\n }\n"
573     if len(Q) > 1: CODE+="}\n"     if len(Q) > 1: CODE+="}\n"
# Line 589  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 589  def makePDE(S, x, Q, W, DIM=2, system=Fa
589          CODE2=""          CODE2=""
590          EM = {}          EM = {}
591          for i in range(len(S)):          for i in range(len(S)):
592          for j in range(len(S)):              for j in range(len(S)):
593              EM[(i,j)] = 0                  EM[(i,j)] = 0
594                                
595          for q in xrange(len(Q)):          for q in range(len(Q)):
596          for dj in range(DIM):              for dj in range(DIM):
597                A_name="C_%d_%d"%(dj,q)                      A_name="C_%d_%d"%(dj,q)
598                    A=Symbol(A_name)                    A=Symbol(A_name)
599                    DATA_C.append(A)                    DATA_C.append(A)
600                    if system:                    if system:
601                CODE2+="const register double %s = C_p[INDEX4(k,m,%s, %s, p.numEqu,p.numComp,%s)];\n"%(A_name, dj,  q, DIM)                        CODE2+="const register double %s = C_p[INDEX4(k,m,%s, %s, p.numEqu,p.numComp,%s)];\n"%(A_name, dj,  q, DIM)
602            else:                    else:
603                CODE2+="const register double %s = C_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj, q, DIM)                        CODE2+="const register double %s = C_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj, q, DIM)
604            for i in range(len(S)):                     for i in range(len(S)):
605                    for j in range(len(S)):                        for j in range(len(S)):
606                        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) ] )                            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) ] )
607          CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system)          CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system)
608                    
609      if system: CODE+="}\n }\n"          if system: CODE+="}\n }\n"
610      CODE+="} else { /* constant data */\n"          CODE+="} else { /* constant data */\n"
611     if system:     if system:
612         CODE+= """for (k=0;k<p.numEqu;k++) {         CODE+= """for (k=0;k<p.numEqu;k++) {
613                              for (m=0;m<p.numComp;m++) {                              for (m=0;m<p.numComp;m++) {
# Line 617  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 617  def makePDE(S, x, Q, W, DIM=2, system=Fa
617     EM = {}     EM = {}
618     for i in range(len(S)):     for i in range(len(S)):
619         for j in range(len(S)):         for j in range(len(S)):
620        EM[(i,j)] = 0            EM[(i,j)] = 0
621     for dj in range(DIM):     for dj in range(DIM):
622           A_name="C_%d"%(dj)                 A_name="C_%d"%(dj)
623               A=Symbol(A_name)               A=Symbol(A_name)
624               DATA_C.append(A)               DATA_C.append(A)
625               if system:               if system:
626                CODE2+="const register double %s = C_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, dj)                        CODE2+="const register double %s = C_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, dj)
627           else:               else:
628                CODE2+="const register double %s = C_p[%s];\n"%(A_name, dj)                        CODE2+="const register double %s = C_p[%s];\n"%(A_name, dj)
629               for q in xrange(len(Q)):               for q in range(len(Q)):
630            for i in range(len(S)):                     for i in range(len(S)):
631                    for j in range(len(S)):                        for j in range(len(S)):
632                        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) ] )                            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) ] )
633     CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system)         CODE+=CODE2+generatePDECode(DATA_C, EM, GLOBAL_TMP,system)    
634     if system:  CODE+="}\n }\n"     if system:  CODE+="}\n }\n"
635     if len(Q) > 1: CODE+="}\n"     if len(Q) > 1: CODE+="}\n"
# Line 651  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 651  def makePDE(S, x, Q, W, DIM=2, system=Fa
651          CODE2=""          CODE2=""
652          EM = {}          EM = {}
653          for i in range(len(S)):          for i in range(len(S)):
654          for j in range(len(S)):              for j in range(len(S)):
655              EM[(i,j)] = 0                  EM[(i,j)] = 0
656                                
657          for q in xrange(len(Q)):          for q in range(len(Q)):
658                A_name="D_%d"%(q,)                      A_name="D_%d"%(q,)
659                    A=Symbol(A_name)                    A=Symbol(A_name)
660                    DATA_D.append(A)                    DATA_D.append(A)
661                    if system:                    if system:
662                CODE2+="const register double %s = D_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, q)                        CODE2+="const register double %s = D_p[INDEX3(k,m,%s, p.numEqu,p.numComp)];\n"%(A_name, q)
663            else:                    else:
664                CODE2+="const register double %s = D_p[%s];\n"%(A_name, q)                        CODE2+="const register double %s = D_p[%s];\n"%(A_name, q)
665            for i in range(len(S)):                     for i in range(len(S)):
666                    for j in range(len(S)):                        for j in range(len(S)):
667                        EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )                            EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] )
668          CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system)          CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system)
669                    
670      if system: CODE+="}\n }\n"          if system: CODE+="}\n }\n"
671      CODE+="} else { /* constant data */\n"          CODE+="} else { /* constant data */\n"
672     if system:     if system:
673         CODE+= """for (k=0;k<p.numEqu;k++) {         CODE+= """for (k=0;k<p.numEqu;k++) {
674                              for (m=0;m<p.numComp;m++) {                              for (m=0;m<p.numComp;m++) {
# Line 678  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 678  def makePDE(S, x, Q, W, DIM=2, system=Fa
678     EM = {}     EM = {}
679     for i in range(len(S)):     for i in range(len(S)):
680         for j in range(len(S)):         for j in range(len(S)):
681        EM[(i,j)] = 0            EM[(i,j)] = 0
682     if 1:     if 1:
683           A_name="D_0"                 A_name="D_0"
684               A=Symbol(A_name)               A=Symbol(A_name)
685               DATA_D.append(A)               DATA_D.append(A)
686               if system:               if system:
687                CODE2+="const register double %s = D_p[INDEX2(k,m, p.numEqu)];\n"%(A_name,)                        CODE2+="const register double %s = D_p[INDEX2(k,m, p.numEqu)];\n"%(A_name,)
688           else:               else:
689                CODE2+="const register double %s = D_p[0];\n"%(A_name,)                        CODE2+="const register double %s = D_p[0];\n"%(A_name,)
690               for q in xrange(len(Q)):               for q in range(len(Q)):
691            for i in range(len(S)):                     for i in range(len(S)):
692                    for j in range(len(S)):                        for j in range(len(S)):
693                        EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i] ).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )                            EM[(i,j)] = EM[(i,j)] + (A * W[q] * S[j] * S[i] ).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] )
694     CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system)         CODE+=CODE2+generatePDECode(DATA_D, EM, GLOBAL_TMP,system)    
695     if system:  CODE+="}\n }\n"     if system:  CODE+="}\n }\n"
696     if len(Q) > 1: CODE+="}\n"     if len(Q) > 1: CODE+="}\n"
# Line 712  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 712  def makePDE(S, x, Q, W, DIM=2, system=Fa
712          CODE2=""          CODE2=""
713          EM = {}          EM = {}
714          for i in range(len(S)):          for i in range(len(S)):
715              EM[i] = 0                  EM[i] = 0
716                                
717          for q in xrange(len(Q)):          for q in range(len(Q)):
718          for dj in range(DIM):              for dj in range(DIM):
719                A_name="X_%d_%d"%(dj,q)                      A_name="X_%d_%d"%(dj,q)
720                    A=Symbol(A_name)                    A=Symbol(A_name)
721                    DATA_X.append(A)                    DATA_X.append(A)
722                    if system:                    if system:
723                CODE2+="const register double %s = X_p[INDEX3(k,%s, %s, p.numEqu,%s)];\n"%(A_name, dj,  q, DIM)                        CODE2+="const register double %s = X_p[INDEX3(k,%s, %s, p.numEqu,%s)];\n"%(A_name, dj,  q, DIM)
724            else:                    else:
725                CODE2+="const register double %s = X_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj,q,DIM)                        CODE2+="const register double %s = X_p[INDEX2(%s,%s,%s)];\n"%(A_name, dj,q,DIM)
726            for j in range(len(S)):                     for j in range(len(S)):
727                        EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )                            EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] )
728          CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system)          CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system)
729                    
730      if system: CODE+="}\n"          if system: CODE+="}\n"
731      CODE+="} else { /* constant data */\n"          CODE+="} else { /* constant data */\n"
732     if system:     if system:
733         CODE+= """for (k=0;k<p.numEqu;k++) {         CODE+= """for (k=0;k<p.numEqu;k++) {
734              """              """
# Line 736  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 736  def makePDE(S, x, Q, W, DIM=2, system=Fa
736     CODE2=""     CODE2=""
737     EM = {}     EM = {}
738     for i in range(len(S)):     for i in range(len(S)):
739        EM[i] = 0            EM[i] = 0
740     for dj in range(DIM):     for dj in range(DIM):
741           A_name="X_%d"%(dj)                 A_name="X_%d"%(dj)
742               A=Symbol(A_name)               A=Symbol(A_name)
743               DATA_X.append(A)               DATA_X.append(A)
744               if system:               if system:
745                CODE2+="const register double %s = X_p[INDEX2(k,%s, p.numEqu)];\n"%(A_name, dj)                        CODE2+="const register double %s = X_p[INDEX2(k,%s, p.numEqu)];\n"%(A_name, dj)
746           else:               else:
747                CODE2+="const register double %s = X_p[%s];\n"%(A_name, dj)                        CODE2+="const register double %s = X_p[%s];\n"%(A_name, dj)
748               for q in xrange(len(Q)):               for q in range(len(Q)):
749            for j in range(len(S)):                     for j in range(len(S)):
750                        EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )                            EM[j] = EM[j] + (A * W[q] * diff(S[j],x[dj])).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] )
751     CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system)         CODE+=CODE2+generatePDECode(DATA_X, EM, GLOBAL_TMP,system)    
752     if system:  CODE+="}\n"     if system:  CODE+="}\n"
753     if len(Q) > 1: CODE+="}\n"     if len(Q) > 1: CODE+="}\n"
# Line 769  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 769  def makePDE(S, x, Q, W, DIM=2, system=Fa
769          CODE2=""          CODE2=""
770          EM = {}          EM = {}
771          for i in range(len(S)):          for i in range(len(S)):
772              EM[i] = 0                  EM[i] = 0
773                                
774          for q in xrange(len(Q)):          for q in range(len(Q)):
775                A_name="Y_%d"%(q,)                      A_name="Y_%d"%(q,)
776                    A=Symbol(A_name)                    A=Symbol(A_name)
777                    DATA_Y.append(A)                    DATA_Y.append(A)
778                    if system:                    if system:
779                CODE2+="const register double %s = Y_p[INDEX3(k,%s, p.numEqu)];\n"%(A_name, q)                        CODE2+="const register double %s = Y_p[INDEX3(k,%s, p.numEqu)];\n"%(A_name, q)
780            else:                    else:
781                CODE2+="const register double %s = Y_p[%s];\n"%(A_name, q)                        CODE2+="const register double %s = Y_p[%s];\n"%(A_name, q)
782            for i in range(len(S)):                     for i in range(len(S)):
783                        EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )                            EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] )
784          CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system)          CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system)
785                    
786      if system: CODE+="}\n"          if system: CODE+="}\n"
787      CODE+="} else { /* constant data */\n"          CODE+="} else { /* constant data */\n"
788     if system:     if system:
789         CODE+= """for (k=0;k<p.numEqu;k++) {         CODE+= """for (k=0;k<p.numEqu;k++) {
790              """              """
# Line 792  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 792  def makePDE(S, x, Q, W, DIM=2, system=Fa
792     CODE2=""     CODE2=""
793     EM = {}     EM = {}
794     for i in range(len(S)):     for i in range(len(S)):
795        EM[i] = 0            EM[i] = 0
796     if 1:     if 1:
797           A_name="Y_0"                 A_name="Y_0"
798               A=Symbol(A_name)               A=Symbol(A_name)
799               DATA_Y.append(A)               DATA_Y.append(A)
800               if system:               if system:
801                CODE2+="const register double %s = Y_p[k];\n"%(A_name,)                        CODE2+="const register double %s = Y_p[k];\n"%(A_name,)
802           else:               else:
803                CODE2+="const register double %s = Y_p[0];\n"%(A_name,)                        CODE2+="const register double %s = Y_p[0];\n"%(A_name,)
804               for q in xrange(len(Q)):               for q in range(len(Q)):
805            for i in range(len(S)):                     for i in range(len(S)):
806                        EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in xrange(DIM) ] )                            EM[i] = EM[i] + (A * W[q] * S[i]).subs( [ (x[jj], Q[q][jj]) for jj in range(DIM) ] )
807     CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system)         CODE+=CODE2+generatePDECode(DATA_Y, EM, GLOBAL_TMP,system)    
808     if system:  CODE+="}\n"     if system:  CODE+="}\n"
809     if len(Q) > 1: CODE+="}\n"     if len(Q) > 1: CODE+="}\n"
Legend:
Removed from v.4575  
changed lines
  Added in v.4576
  ViewVC Help
Powered by ViewVC 1.1.26