/[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 3791 by gross, Wed Nov 23 01:39:45 2011 UTC revision 3792 by caltinay, Wed Feb 1 06:16:25 2012 UTC
# Line 10  DIGITS=20 Line 10  DIGITS=20
10  #  #
11  #  quadrature form in [0,1]  #  quadrature form in [0,1]
12  #  #
13  q0=sqrt(RealNumber(15))/5  
14  qD1=[ (1-q0)/2,  RealNumber(1)/2,  (1+q0)/2 ]  q0=1/sqrt(RealNumber(3))
15  wD1=[RealNumber(5)/18, RealNumber(4)/ 9, RealNumber(5)/18 ]  qD1=[ (1-q0)/2,  (1+q0)/2 ]
16    wD1=[RealNumber(1)/2, RealNumber(1)/2 ]
17    
18    
19  q0=1/sqrt(RealNumber(3))  q0=1/sqrt(RealNumber(3))
20  qD1=[ (1-q0)/2,  (1+q0)/2 ]  qD1=[ (1-q0)/2,  (1+q0)/2 ]
# Line 127  def generate(DIM): Line 129  def generate(DIM):
129    
130     # interpolation to Quad points     # interpolation to Quad points
131     CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"     CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"
    CODE+=createIntegrationCode(Q, W, loopindex=idx_full)  
    CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"  
    CODE+=createIntegrationCode(Q_r, W_r, loopindex=idx_full)  
    CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n"  
    for i in xrange(len(Q_r_faces)):  
         CODE+="if (face_offset(%s)>-1) {\n"%i  
         CODE+=createIntegrationCode(Q_faces[i], W_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)  
         CODE+="\n} /* end of face %s */\n"%i  
    CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n"  
    for i in xrange(len(Q_r_faces)):  
         CODE+="if (face_offset(%s)>-1) {\n"%i  
         CODE+=createIntegrationCode(Q_r_faces[i], W_r_faces[i], loopindex=idx_faces[i],  gridoffset="face_offset(%s)"%i)  
         CODE+="\n} /* end of face %s */\n"%i  
    CODE+="\n} /* end of out_data_type branching */\n"  
    insertCode("Assemble_Integration_%sD.c"%DIM, { "SNIP" : CODE})  
    1/0  
   
    # interpolation to Quad points  
    CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"  
132     CODE+=createCode(f_interpolation, x, Q, loopindex=idx_full)     CODE+=createCode(f_interpolation, x, Q, loopindex=idx_full)
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)
# Line 161  def generate(DIM): Line 144  def generate(DIM):
144          CODE+="\n} /* end of face %s */\n"%i          CODE+="\n} /* end of face %s */\n"%i
145     CODE+="\n} /* end of out_data_type branching */\n"     CODE+="\n} /* end of out_data_type branching */\n"
146     insertCode("Assemble_Interpolation_%sD.c"%DIM, { "SNIP" : CODE})     insertCode("Assemble_Interpolation_%sD.c"%DIM, { "SNIP" : CODE})
147      
148          # gradient to Quad points
    # gradient to Quad points  
149     CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"     CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n"
150     CODE+=createGradientCode(f_interpolation, x, Q, loopindex=idx_full)     CODE+=createGradientCode(f_interpolation, x, Q, loopindex=idx_full)
151     CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"     CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n"
# Line 180  def generate(DIM): Line 162  def generate(DIM):
162          CODE+="\n} /* end of face %s */\n"%i          CODE+="\n} /* end of face %s */\n"%i
163     CODE+="\n} /* end of out_data_type branching */\n"     CODE+="\n} /* end of out_data_type branching */\n"
164     insertCode("Assemble_Gradient_%sD.c"%DIM, { "SNIP" : CODE})     insertCode("Assemble_Gradient_%sD.c"%DIM, { "SNIP" : CODE})
165      
166     #generate PDE assemblage     #generate PDE assemblage
167     CODE,PRECODE = makePDE(S, x, Q, W, DIM=DIM, system=True)     CODE,PRECODE = makePDE(S, x, Q, W, DIM=DIM, system=True)
168     insertCode("Assemble_PDE_System_%sD.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE } )     insertCode("Assemble_PDE_System_%sD.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE } )
# Line 404  def createCode(F, x, Q, gridoffset="", l Line 386  def createCode(F, x, Q, gridoffset="", l
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        
 def createIntegrationCode(Q, W, gridoffset="", loopindex=(0,1,2)):  
    DIM=len(loopindex)  
    GLOBAL_TMP={}  
    LOCAL_TMP_E={}  
    TXT_E=""  
    for q in range(len(W)):  
         if not GLOBAL_TMP.has_key(W[q]):  
          GLOBAL_TMP[W[q]]=Symbol("w_%s"%len(GLOBAL_TMP))  
         A=Symbol("f_%s"%q)  
         TXT_E+="const register double %s = in[INDEX3(i,%s,e, NCOMP,%s)];\n"%(A,q,len(W))  
         if not LOCAL_TMP_E.has_key(GLOBAL_TMP[W[q]]):  
          LOCAL_TMP_E[GLOBAL_TMP[W[q]]]=0  
     LOCAL_TMP_E[GLOBAL_TMP[W[q]]]=LOCAL_TMP_E[GLOBAL_TMP[W[q]]]+A  
    for p in LOCAL_TMP_E.keys():  
         TXT_E+="const register double %s = %s;\n"%( p, ccode(LOCAL_TMP_E[p]))  
    print TXT_E  
    print GLOBAL_TMP  
    1/0  
    #=================  
   
   
     
    k3=""  
    for i in xrange(DIM-1,-1,-1):  
           if loopindex[i] > -1: k3+=",k%i"%loopindex[i]  
    TXT=""  
    for a,v in  consts.items():  
       TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS)))  
    TXT+="#pragma omp parallel for private(i%s)"%k3  
    for i in xrange(DIM-1,-1,-1):  
           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])  
    TXT+="\nfor (i =0; i < NCOMP; ++i) {\n"  
    print TXT  
    1/0  
    #  interpolation to quadrature points  
    for q in xrange(len(Q)):  
       IDX2="INDEX%s(%s,%s)"%(DIM,k,N)  
       if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2  
       TXT+="out[INDEX3(i,%s,%s,NCOMP,%s)] = %s;\n"%(q,IDX2,len(Q),ccode(F[q]))  
    TXT+="} /* close component loop i */\n"  
    for i in xrange(DIM):  
            if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i]  
    return TXT  
     
389  def subPoint(g,x, p):  def subPoint(g,x, p):
390       out=g       out=g
391       for i in range(len(x)): out=out.subs(x[i], p[i])       for i in range(len(x)): out=out.subs(x[i], p[i])
# Line 875  def makePDE(S, x, Q, W, DIM=2, system=Fa Line 813  def makePDE(S, x, Q, W, DIM=2, system=Fa
813         PRECODE+="const double %s = %s;\n"%(GLOBAL_TMP[t],ccode(t.evalf(n=DIGITS)))         PRECODE+="const double %s = %s;\n"%(GLOBAL_TMP[t],ccode(t.evalf(n=DIGITS)))
814     return CODE, PRECODE     return CODE, PRECODE
815    
816  for d in [2]:  for d in [3]:
817       generate(d)       generate(d)

Legend:
Removed from v.3791  
changed lines
  Added in v.3792

  ViewVC Help
Powered by ViewVC 1.1.26