/[escript]/trunk/escript/py_src/generatediff
ViewVC logotype

Diff of /trunk/escript/py_src/generatediff

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

revision 438 by gross, Fri Jan 20 00:39:00 2006 UTC revision 441 by gross, Fri Jan 20 03:40:39 2006 UTC
# Line 696  def minimumTEST(arg0,arg1): Line 696  def minimumTEST(arg0,arg1):
696                else:                else:
697                 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]                 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
698       return out       return out
699  def unrollLoops(a,b,o,arg,tap=""):  def unrollLoops(a,b,o,arg,tap="",x="x"):
700      out=""      out=""
701      if a.rank==1:      if a.rank==1:
702               z=""               z=""
703               for i99 in range(a.shape[0]):               for i99 in range(a.shape[0]):
704                 if not z=="": z+="+"                 if not z=="": z+="+"
705                 if o=="1":                 if o=="1":
706                    z+="(%s)*x[%s]"%(a[i99]+b[i99],i99)                    z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
707                 else:                 else:
708                    z+="(%s)*x[%s]**o+(%s)*x[%s]"%(a[i99],i99,b[i99],i99)                    z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
709    
710               out+=tap+"%s=%s\n"%(arg,z)               out+=tap+"%s=%s\n"%(arg,z)
711    
# Line 717  def unrollLoops(a,b,o,arg,tap=""): Line 717  def unrollLoops(a,b,o,arg,tap=""):
717                 if o=="1":                 if o=="1":
718                    z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)                    z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
719                 else:                 else:
720                    z+="(%s)*x[%s]**o+(%s)*x[%s]"%(a[i0,i99],i99,b[i0,i99],i99)                    z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
721    
722               out+=tap+"%s[%s]=%s\n"%(arg,i0,z)               out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
723      elif a.rank==3:      elif a.rank==3:
# Line 727  def unrollLoops(a,b,o,arg,tap=""): Line 727  def unrollLoops(a,b,o,arg,tap=""):
727               for i99 in range(a.shape[2]):               for i99 in range(a.shape[2]):
728                 if not z=="": z+="+"                 if not z=="": z+="+"
729                 if o=="1":                 if o=="1":
730                    z+="(%s)*x[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],i99)                    z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
731                 else:                 else:
732                    z+="(%s)*x[%s]**o+(%s)*x[%s]"%(a[i0,i1,i99],i99,b[i0,i1,i99],i99)                    z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
733    
734               out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)               out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
735      elif a.rank==4:      elif a.rank==4:
# Line 740  def unrollLoops(a,b,o,arg,tap=""): Line 740  def unrollLoops(a,b,o,arg,tap=""):
740               for i99 in range(a.shape[3]):               for i99 in range(a.shape[3]):
741                 if not z=="": z+="+"                 if not z=="": z+="+"
742                 if o=="1":                 if o=="1":
743                    z+="(%s)*x[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],i99)                    z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
744                 else:                 else:
745                    z+="(%s)*x[%s]**o+(%s)*x[%s]"%(a[i0,i1,i2,i99],i99,b[i0,i1,i2,i99],i99)                    z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
746    
747               out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)               out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
748      elif a.rank==5:      elif a.rank==5:
# Line 754  def unrollLoops(a,b,o,arg,tap=""): Line 754  def unrollLoops(a,b,o,arg,tap=""):
754               for i99 in range(a.shape[4]):               for i99 in range(a.shape[4]):
755                 if not z=="": z+="+"                 if not z=="": z+="+"
756                 if o=="1":                 if o=="1":
757                    z+="(%s)*x[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],i99)                    z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
758                 else:                 else:
759                    z+="(%s)*x[%s]**o+(%s)*x[%s]"%(a[i0,i1,i2,i3,i99],i99,b[i0,i1,i2,i3,i99],i99)                    z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
760    
761               out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)               out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
762      return out      return out
# Line 797  def unrollLoopsOfGrad(a,b,o,arg,tap=""): Line 797  def unrollLoopsOfGrad(a,b,o,arg,tap=""):
797                 else:                 else:
798                   out+=tap+"%s[%s,%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99],i99,b[i0,i1,i2,i99])                   out+=tap+"%s[%s,%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99],i99,b[i0,i1,i2,i99])
799      return out      return out
800    def unrollLoopsOfDiv(a,b,o,arg,tap=""):
801    
802    
803        
804        out=tap+arg+"="
805        if o=="1":
806           z=0.
807           for i99 in range(a.shape[0]):
808                z+=b[i99,i99]+a[i99,i99]
809           out+="(%s)"%z    
810        else:
811           z=0.
812           for i99 in range(a.shape[0]):
813                z+=b[i99,i99]
814                if i99>0: out+="+"
815                out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
816           out+="+(%s)"%z    
817        return out
818    
819  def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):  def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
820      if where=="Function":      if where=="Function":
821         xfac_o=1.         xfac_o=1.
# Line 916  def unrollLoopsOfInteriorIntegral(a,b,wh Line 935  def unrollLoopsOfInteriorIntegral(a,b,wh
935                 out+="+(%s)*0.5**o\n"%zop                 out+="+(%s)*0.5**o\n"%zop
936    
937      return out      return out
938    #=======================================================================================================
939    # div
940    #=======================================================================================================
941    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
942      for data in ["Data","Symbol"]:
943         for case in ["ContinuousFunction","Solution","ReducedSolution"]:
944             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
945             tname="test_div_on%s_from%s_%s"%(where,data,case)
946             text+="   def %s(self):\n"%tname
947             text+="      \"\"\"\n"
948             text+="      tests divergence of %s on the %s\n\n"%(data,where)
949             text+="      assumptions: %s(self.domain) exists\n"%case
950             text+="                   self.domain supports div on %s\n"%where
951             text+="      \"\"\"\n"
952             if case=="ReducedSolution":
953                text+="      o=1\n"
954                o="1"
955             else:
956                text+="      o=self.order\n"
957                o="o"
958             text+="      dim=self.domain.getDim()\n"
959             text+="      w_ref=%s(self.domain)\n"%where
960             text+="      x_ref=w_ref.getX()\n"
961             text+="      w=%s(self.domain)\n"%case
962             text+="      x=w.getX()\n"
963             a_2=makeArray((2,2),[-1.,1])
964             b_2=makeArray((2,2),[-1.,1])
965             a_3=makeArray((3,3),[-1.,1])
966             b_3=makeArray((3,3),[-1.,1])
967             if data=="Symbol":
968                text+="      arg=Symbol(shape=(dim,),dim=dim)\n"
969                val="s"
970                res="sub"
971             else:
972                val="arg"
973                res="res"
974             text+="      %s=Vector(0,w)\n"%val
975             text+="      if dim==2:\n"
976             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
977             text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap="        ")
978             text+="\n      else:\n"
979            
980             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
981             text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap="        ")
982             text+="\n      res=div(arg,where=w_ref)\n"
983             if data=="Symbol":
984                text+="      sub=res.substitute({arg:s})\n"
985             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
986             text+="      self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
987             text+="      self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
988             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
989    
990    
991             t_prog+=text
992    print t_prog
993    1/0
994    
995    #=======================================================================================================
996    # interpolation
997    #=======================================================================================================
998    for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
999      for data in ["Data","Symbol"]:
1000         for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1001          for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1002            if  where==case or \
1003                ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
1004                ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
1005                (case=="ContinuousFunction" and  where in ["Solution","ReducedSolution"]) or \
1006                (case=="Solution" and  where=="ReducedSolution") :
1007                
1008    
1009             text="   #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1010             tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1011             text+="   def %s(self):\n"%tname
1012             text+="      \"\"\"\n"
1013             text+="      tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
1014             text+="      assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
1015             text+="      \"\"\"\n"
1016             if case=="ReducedSolution" or where=="ReducedSolution":
1017                text+="      o=1\n"
1018                o="1"
1019             else:
1020                text+="      o=self.order\n"
1021                o="o"
1022             text+="      dim=self.domain.getDim()\n"
1023             text+="      w_ref=%s(self.domain)\n"%where
1024             text+="      x_ref=w_ref.getX()\n"
1025             text+="      w=%s(self.domain)\n"%case
1026             text+="      x=w.getX()\n"
1027             a_2=makeArray(sh+(2,),[-1.,1])
1028             b_2=makeArray(sh+(2,),[-1.,1])
1029             a_3=makeArray(sh+(3,),[-1.,1])
1030             b_3=makeArray(sh+(3,),[-1.,1])
1031             if data=="Symbol":
1032                text+="      arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1033                val="s"
1034                res="sub"
1035             else:
1036                val="arg"
1037                res="res"
1038             text+="      %s=Data(0,%s,w)\n"%(val,str(sh))
1039             text+="      ref=Data(0,%s,w_ref)\n"%str(sh)
1040             text+="      if dim==2:\n"
1041             text+=unrollLoops(a_2,b_2,o,val,tap="        ")
1042             text+=unrollLoops(a_2,b_2,o,"ref",tap="        ",x="x_ref")
1043             text+="      else:\n"
1044            
1045             text+=unrollLoops(a_3,b_3,o,val,tap="        ")
1046             text+=unrollLoops(a_3,b_3,o,"ref",tap="        ",x="x_ref")
1047             text+="      res=interpolate(arg,where=w_ref)\n"
1048             if data=="Symbol":
1049                text+="      sub=res.substitute({arg:s})\n"
1050             text+="      self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1051             text+="      self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
1052             text+="      self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1053             text+="      self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1054             t_prog+=text
1055    print test_header
1056    print t_prog
1057    print test_tail          
1058    1/0
1059    
1060  #=======================================================================================================  #=======================================================================================================
1061  # grad  # grad

Legend:
Removed from v.438  
changed lines
  Added in v.441

  ViewVC Help
Powered by ViewVC 1.1.26