/[escript]/trunk/escript/py_src/util.py
ViewVC logotype

Diff of /trunk/escript/py_src/util.py

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

revision 396 by gross, Wed Dec 21 05:08:25 2005 UTC revision 536 by gross, Fri Feb 17 03:20:53 2006 UTC
# Line 30  __date__="$Date$" Line 30  __date__="$Date$"
30    
31  import math  import math
32  import numarray  import numarray
33    import numarray.linear_algebra
34  import escript  import escript
35  import os  import os
36    
# Line 43  import os Line 44  import os
44  # def matchType(arg0=0.,arg1=0.):  # def matchType(arg0=0.,arg1=0.):
45  # def matchShape(arg0,arg1):  # def matchShape(arg0,arg1):
46    
 # def maximum(arg0,arg1):  
 # def minimum(arg0,arg1):  
 # def clip(arg,minval,maxval)  
   
 # def transpose(arg,axis=None):  
 # def trace(arg,axis0=0,axis1=1):  
47  # def reorderComponents(arg,index):  # def reorderComponents(arg,index):
48    
 # def integrate(arg,where=None):  
 # def interpolate(arg,where):  
 # def div(arg,where=None):  
 # def grad(arg,where=None):  
   
49  #  #
50  # slicing: get  # slicing: get
51  #          set  #          set
# Line 126  def kronecker(d=3): Line 116  def kronecker(d=3):
116     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
117    
118     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
119     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
120     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
121     @rtype d: L{numarray.NumArray} of rank 2.     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2.
    @remark: the function is identical L{identity}  
122     """     """
123     return identityTensor(d)     return identityTensor(d)
124    
# Line 148  def identity(shape=()): Line 137  def identity(shape=()):
137        if len(shape)==1:        if len(shape)==1:
138            for i0 in range(shape[0]):            for i0 in range(shape[0]):
139               out[i0,i0]=1.               out[i0,i0]=1.
   
140        elif len(shape)==2:        elif len(shape)==2:
141            for i0 in range(shape[0]):            for i0 in range(shape[0]):
142               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 164  def identityTensor(d=3): Line 152  def identityTensor(d=3):
152     return the dxd identity matrix     return the dxd identity matrix
153    
154     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
155     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
156     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
157     @rtype: L{numarray.NumArray} of rank 2.     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2
158     """     """
159     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
160        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
161     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
162           return identity((d.getDim(),))
163       else:
164           return identity((d,))
165    
166  def identityTensor4(d=3):  def identityTensor4(d=3):
167     """     """
# Line 179  def identityTensor4(d=3): Line 170  def identityTensor4(d=3):
170     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
171     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
172     @return: the object u of rank 4 with M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise     @return: the object u of rank 4 with M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise
173     @rtype: L{numarray.NumArray} of rank 4.     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 4.
174     """     """
175     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
176        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
177     return identity((d,d))     elif isinstance(d,escript.Domain):
178           return identity((d.getDim(),d.getDim()))
179       else:
180           return identity((d,d))
181    
182  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
183     """     """
# Line 192  def unitVector(i=0,d=3): Line 186  def unitVector(i=0,d=3):
186     @param i: index     @param i: index
187     @type i: C{int}     @type i: C{int}
188     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
189     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
190     @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise     @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise
191     @rtype: L{numarray.NumArray} of rank 1.     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 1
192     """     """
193     return kronecker(d)[i]     return kronecker(d)[i]
194    
# Line 834  class Symbol(object): Line 828  class Symbol(object):
828         """         """
829         return power(other,self)         return power(other,self)
830    
831       def __getitem__(self,index):
832           """
833           returns the slice defined by index
834    
835           @param index: defines a
836           @type index: C{slice} or C{int} or a C{tuple} of them
837           @return: a S{Symbol} representing the slice defined by index
838           @rtype: L{DependendSymbol}
839           """
840           return GetSlice_Symbol(self,index)
841    
842  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
843     """     """
844     DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.     DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.
# Line 884  class DependendSymbol(Symbol): Line 889  class DependendSymbol(Symbol):
889  #=========================================================  #=========================================================
890  #  Unary operations prserving the shape  #  Unary operations prserving the shape
891  #========================================================  #========================================================
892    class GetSlice_Symbol(DependendSymbol):
893       """
894       L{Symbol} representing getting a slice for a L{Symbol}
895       """
896       def __init__(self,arg,index):
897          """
898          initialization of wherePositive L{Symbol} with argument arg
899          @param arg: argument
900          @type arg: L{Symbol}.
901          @param index: defines index
902          @type index: C{slice} or C{int} or a C{tuple} of them
903          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
904          @raises ValueError: if a step is given
905          """
906          if not isinstance(index,tuple): index=(index,)
907          if len(index)>arg.getRank():
908               raise IndexError,"GetSlice_Symbol: index out of range."
909          sh=()
910          index2=()
911          for i in range(len(index)):
912             ix=index[i]
913             if isinstance(ix,int):
914                if ix<0 or ix>=arg.getShape()[i]:
915                   raise ValueError,"GetSlice_Symbol: index out of range."
916                index2=index2+(ix,)
917             else:
918               if not ix.step==None:
919                 raise ValueError,"GetSlice_Symbol: steping is not supported."
920               if ix.start==None:
921                  s=0
922               else:
923                  s=ix.start
924               if ix.stop==None:
925                  e=arg.getShape()[i]
926               else:
927                  e=ix.stop
928                  if e>arg.getShape()[i]:
929                     raise IndexError,"GetSlice_Symbol: index out of range."
930               index2=index2+(slice(s,e),)
931               if e>s:
932                   sh=sh+(e-s,)
933               elif s>e:
934                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
935          for i in range(len(index),arg.getRank()):
936              index2=index2+(slice(0,arg.getShape()[i]),)
937              sh=sh+(arg.getShape()[i],)
938          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
939    
940       def getMyCode(self,argstrs,format="escript"):
941          """
942          returns a program code that can be used to evaluate the symbol.
943    
944          @param argstrs: gives for each argument a string representing the argument for the evaluation.
945          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
946          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
947          @type format: C{str}
948          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
949          @rtype: C{str}
950          @raise: NotImplementedError: if the requested format is not available
951          """
952          if format=="escript" or format=="str"  or format=="text":
953             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
954          else:
955             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
956    
957       def substitute(self,argvals):
958          """
959          assigns new values to symbols in the definition of the symbol.
960          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
961    
962          @param argvals: new values assigned to symbols
963          @type argvals: C{dict} with keywords of type L{Symbol}.
964          @return: result of the substitution process. Operations are executed as much as possible.
965          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
966          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
967          """
968          if argvals.has_key(self):
969             arg=argvals[self]
970             if self.isAppropriateValue(arg):
971                return arg
972             else:
973                raise TypeError,"%s: new value is not appropriate."%str(self)
974          else:
975             args=self.getSubstitutedArguments(argvals)
976             arg=args[0]
977             index=args[1]
978             return arg.__getitem__(index)
979    
980  def log10(arg):  def log10(arg):
981     """     """
982     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
# Line 916  def wherePositive(arg): Line 1009  def wherePositive(arg):
1009     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1010     """     """
1011     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1012        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))*1.
1013        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out)
1014        return out        return out
1015     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
# Line 998  def whereNegative(arg): Line 1091  def whereNegative(arg):
1091     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1092     """     """
1093     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1094        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))*1.
1095        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out)
1096        return out        return out
1097     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
# Line 1080  def whereNonNegative(arg): Line 1173  def whereNonNegative(arg):
1173     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1174     """     """
1175     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1176        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.
1177        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out)
1178        return out        return out
1179     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
# Line 2862  def length(arg): Line 2955  def length(arg):
2955     """     """
2956     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2957    
2958    def trace(arg,axis_offset=0):
2959       """
2960       returns the trace of arg which the sum of arg[k,k] over k.
2961    
2962       @param arg: argument
2963       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2964       @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
2965                      axis_offset and axis_offset+1 must be equal.
2966       @type axis_offset: C{int}
2967       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2968       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2969       """
2970       if isinstance(arg,numarray.NumArray):
2971          sh=arg.shape
2972          if len(sh)<2:
2973            raise ValueError,"trace: rank of argument must be greater than 1"
2974          if axis_offset<0 or axis_offset>len(sh)-2:
2975            raise ValueError,"trace: axis_offset must be between 0 and %s"%len(sh)-2
2976          s1=1
2977          for i in range(axis_offset): s1*=sh[i]
2978          s2=1
2979          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2980          if not sh[axis_offset] == sh[axis_offset+1]:
2981            raise ValueError,"trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2982          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
2983          out=numarray.zeros([s1,s2],numarray.Float)
2984          for i1 in range(s1):
2985            for i2 in range(s2):
2986                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
2987          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2988          return out
2989       elif isinstance(arg,escript.Data):
2990          return escript_trace(arg,axis_offset)
2991       elif isinstance(arg,float):
2992          raise TypeError,"trace: illegal argument type float."
2993       elif isinstance(arg,int):
2994          raise TypeError,"trace: illegal argument type int."
2995       elif isinstance(arg,Symbol):
2996          return Trace_Symbol(arg,axis_offset)
2997       else:
2998          raise TypeError,"trace: Unknown argument type."
2999    
3000    def escript_trace(arg,axis_offset): # this should be escript._trace
3001          "arg si a Data objects!!!"
3002          if arg.getRank()<2:
3003            raise ValueError,"escript_trace: rank of argument must be greater than 1"
3004          if axis_offset<0 or axis_offset>arg.getRank()-2:
3005            raise ValueError,"escript_trace: axis_offset must be between 0 and %s"%arg.getRank()-2
3006          s=list(arg.getShape())        
3007          if not s[axis_offset] == s[axis_offset+1]:
3008            raise ValueError,"escript_trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3009          out=escript.Data(0.,tuple(s[0:axis_offset]+s[axis_offset+2:]),arg.getFunctionSpace())
3010          if arg.getRank()==2:
3011             for i0 in range(s[0]):
3012                out+=arg[i0,i0]
3013          elif arg.getRank()==3:
3014             if axis_offset==0:
3015                for i0 in range(s[0]):
3016                      for i2 in range(s[2]):
3017                             out[i2]+=arg[i0,i0,i2]
3018             elif axis_offset==1:
3019                for i0 in range(s[0]):
3020                   for i1 in range(s[1]):
3021                             out[i0]+=arg[i0,i1,i1]
3022          elif arg.getRank()==4:
3023             if axis_offset==0:
3024                for i0 in range(s[0]):
3025                      for i2 in range(s[2]):
3026                         for i3 in range(s[3]):
3027                             out[i2,i3]+=arg[i0,i0,i2,i3]
3028             elif axis_offset==1:
3029                for i0 in range(s[0]):
3030                   for i1 in range(s[1]):
3031                         for i3 in range(s[3]):
3032                             out[i0,i3]+=arg[i0,i1,i1,i3]
3033             elif axis_offset==2:
3034                for i0 in range(s[0]):
3035                   for i1 in range(s[1]):
3036                      for i2 in range(s[2]):
3037                             out[i0,i1]+=arg[i0,i1,i2,i2]
3038          return out
3039    class Trace_Symbol(DependendSymbol):
3040       """
3041       L{Symbol} representing the result of the trace function
3042       """
3043       def __init__(self,arg,axis_offset=0):
3044          """
3045          initialization of trace L{Symbol} with argument arg
3046          @param arg: argument of function
3047          @type arg: L{Symbol}.
3048          @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
3049                      axis_offset and axis_offset+1 must be equal.
3050          @type axis_offset: C{int}
3051          """
3052          if arg.getRank()<2:
3053            raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"
3054          if axis_offset<0 or axis_offset>arg.getRank()-2:
3055            raise ValueError,"Trace_Symbol: axis_offset must be between 0 and %s"%arg.getRank()-2
3056          s=list(arg.getShape())        
3057          if not s[axis_offset] == s[axis_offset+1]:
3058            raise ValueError,"Trace_Symbol: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3059          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3060    
3061       def getMyCode(self,argstrs,format="escript"):
3062          """
3063          returns a program code that can be used to evaluate the symbol.
3064    
3065          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3066          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3067          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3068          @type format: C{str}
3069          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3070          @rtype: C{str}
3071          @raise: NotImplementedError: if the requested format is not available
3072          """
3073          if format=="escript" or format=="str"  or format=="text":
3074             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3075          else:
3076             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3077    
3078       def substitute(self,argvals):
3079          """
3080          assigns new values to symbols in the definition of the symbol.
3081          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3082    
3083          @param argvals: new values assigned to symbols
3084          @type argvals: C{dict} with keywords of type L{Symbol}.
3085          @return: result of the substitution process. Operations are executed as much as possible.
3086          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3087          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3088          """
3089          if argvals.has_key(self):
3090             arg=argvals[self]
3091             if self.isAppropriateValue(arg):
3092                return arg
3093             else:
3094                raise TypeError,"%s: new value is not appropriate."%str(self)
3095          else:
3096             arg=self.getSubstitutedArguments(argvals)
3097             return trace(arg[0],axis_offset=arg[1])
3098    
3099       def diff(self,arg):
3100          """
3101          differential of this object
3102    
3103          @param arg: the derivative is calculated with respect to arg
3104          @type arg: L{escript.Symbol}
3105          @return: derivative with respect to C{arg}
3106          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3107          """
3108          if arg==self:
3109             return identity(self.getShape())
3110          else:
3111             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3112    
3113    def transpose(arg,axis_offset=None):
3114       """
3115       returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.
3116    
3117       @param arg: argument
3118       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3119       @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3120                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3121       @type axis_offset: C{int}
3122       @return: transpose of arg
3123       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3124       """
3125       if isinstance(arg,numarray.NumArray):
3126          if axis_offset==None: axis_offset=int(arg.rank/2)
3127          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3128       elif isinstance(arg,escript.Data):
3129          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3130          return escript_transpose(arg,axis_offset)
3131       elif isinstance(arg,float):
3132          if not ( axis_offset==0 or axis_offset==None):
3133            raise ValueError,"transpose: axis_offset must be 0 for float argument"
3134          return arg
3135       elif isinstance(arg,int):
3136          if not ( axis_offset==0 or axis_offset==None):
3137            raise ValueError,"transpose: axis_offset must be 0 for int argument"
3138          return float(arg)
3139       elif isinstance(arg,Symbol):
3140          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3141          return Transpose_Symbol(arg,axis_offset)
3142       else:
3143          raise TypeError,"transpose: Unknown argument type."
3144    
3145    def escript_transpose(arg,axis_offset): # this should be escript._transpose
3146          "arg si a Data objects!!!"
3147          r=arg.getRank()
3148          if axis_offset<0 or axis_offset>r:
3149            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3150          s=arg.getShape()
3151          s_out=s[axis_offset:]+s[:axis_offset]
3152          out=escript.Data(0.,s_out,arg.getFunctionSpace())
3153          if r==4:
3154             if axis_offset==1:
3155                for i0 in range(s_out[0]):
3156                   for i1 in range(s_out[1]):
3157                      for i2 in range(s_out[2]):
3158                         for i3 in range(s_out[3]):
3159                             out[i0,i1,i2,i3]=arg[i3,i0,i1,i2]
3160             elif axis_offset==2:
3161                for i0 in range(s_out[0]):
3162                   for i1 in range(s_out[1]):
3163                      for i2 in range(s_out[2]):
3164                         for i3 in range(s_out[3]):
3165                             out[i0,i1,i2,i3]=arg[i2,i3,i0,i1]
3166             elif axis_offset==3:
3167                for i0 in range(s_out[0]):
3168                   for i1 in range(s_out[1]):
3169                      for i2 in range(s_out[2]):
3170                         for i3 in range(s_out[3]):
3171                             out[i0,i1,i2,i3]=arg[i1,i2,i3,i0]
3172             else:
3173                for i0 in range(s_out[0]):
3174                   for i1 in range(s_out[1]):
3175                      for i2 in range(s_out[2]):
3176                         for i3 in range(s_out[3]):
3177                             out[i0,i1,i2,i3]=arg[i0,i1,i2,i3]
3178          elif r==3:
3179             if axis_offset==1:
3180                for i0 in range(s_out[0]):
3181                   for i1 in range(s_out[1]):
3182                      for i2 in range(s_out[2]):
3183                             out[i0,i1,i2]=arg[i2,i0,i1]
3184             elif axis_offset==2:
3185                for i0 in range(s_out[0]):
3186                   for i1 in range(s_out[1]):
3187                      for i2 in range(s_out[2]):
3188                             out[i0,i1,i2]=arg[i1,i2,i0]
3189             else:
3190                for i0 in range(s_out[0]):
3191                   for i1 in range(s_out[1]):
3192                      for i2 in range(s_out[2]):
3193                             out[i0,i1,i2]=arg[i0,i1,i2]
3194          elif r==2:
3195             if axis_offset==1:
3196                for i0 in range(s_out[0]):
3197                   for i1 in range(s_out[1]):
3198                             out[i0,i1]=arg[i1,i0]
3199             else:
3200                for i0 in range(s_out[0]):
3201                   for i1 in range(s_out[1]):
3202                             out[i0,i1]=arg[i0,i1]
3203          elif r==1:
3204              for i0 in range(s_out[0]):
3205                   out[i0]=arg[i0]
3206          elif r==0:
3207                 out=arg+0.
3208          return out
3209    class Transpose_Symbol(DependendSymbol):
3210       """
3211       L{Symbol} representing the result of the transpose function
3212       """
3213       def __init__(self,arg,axis_offset=None):
3214          """
3215          initialization of transpose L{Symbol} with argument arg
3216    
3217          @param arg: argument of function
3218          @type arg: L{Symbol}.
3219           @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3220                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3221          @type axis_offset: C{int}
3222          """
3223          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3224          if axis_offset<0 or axis_offset>arg.getRank():
3225            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3226          s=arg.getShape()
3227          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3228    
3229       def getMyCode(self,argstrs,format="escript"):
3230          """
3231          returns a program code that can be used to evaluate the symbol.
3232    
3233          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3234          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3235          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3236          @type format: C{str}
3237          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3238          @rtype: C{str}
3239          @raise: NotImplementedError: if the requested format is not available
3240          """
3241          if format=="escript" or format=="str"  or format=="text":
3242             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3243          else:
3244             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3245    
3246       def substitute(self,argvals):
3247          """
3248          assigns new values to symbols in the definition of the symbol.
3249          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3250    
3251          @param argvals: new values assigned to symbols
3252          @type argvals: C{dict} with keywords of type L{Symbol}.
3253          @return: result of the substitution process. Operations are executed as much as possible.
3254          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3255          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3256          """
3257          if argvals.has_key(self):
3258             arg=argvals[self]
3259             if self.isAppropriateValue(arg):
3260                return arg
3261             else:
3262                raise TypeError,"%s: new value is not appropriate."%str(self)
3263          else:
3264             arg=self.getSubstitutedArguments(argvals)
3265             return transpose(arg[0],axis_offset=arg[1])
3266    
3267       def diff(self,arg):
3268          """
3269          differential of this object
3270    
3271          @param arg: the derivative is calculated with respect to arg
3272          @type arg: L{escript.Symbol}
3273          @return: derivative with respect to C{arg}
3274          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3275          """
3276          if arg==self:
3277             return identity(self.getShape())
3278          else:
3279             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3280    def symmetric(arg):
3281        """
3282        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3283    
3284        @param arg: square matrix. Must have rank 2 or 4 and be square.
3285        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3286        @return: symmetric part of arg
3287        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3288        """
3289        return (arg+transpose(arg))/2
3290    
3291    def nonsymmetric(arg):
3292        """
3293        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3294    
3295        @param arg: square matrix. Must have rank 2 or 4 and be square.
3296        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3297        @return: symmetric part of arg
3298        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3299        """
3300        return (arg-transpose(arg))/2
3301    
3302    def inverse(arg):
3303        """
3304        returns the inverse of the square matrix arg.
3305    
3306        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3307        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3308        @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3309        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3310        @remark: for L{escript.Data} objects the dimension is restricted to 3.
3311        """
3312        if isinstance(arg,numarray.NumArray):
3313          return numarray.linear_algebra.inverse(arg)
3314        elif isinstance(arg,escript.Data):
3315          return escript_inverse(arg)
3316        elif isinstance(arg,float):
3317          return 1./arg
3318        elif isinstance(arg,int):
3319          return 1./float(arg)
3320        elif isinstance(arg,Symbol):
3321          return Inverse_Symbol(arg)
3322        else:
3323          raise TypeError,"inverse: Unknown argument type."
3324    
3325    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3326          "arg is a Data objects!!!"
3327          if not arg.getRank()==2:
3328            raise ValueError,"escript_inverse: argument must have rank 2"
3329          s=arg.getShape()      
3330          if not s[0] == s[1]:
3331            raise ValueError,"escript_inverse: argument must be a square matrix."
3332          out=escript.Data(0.,s,arg.getFunctionSpace())
3333          if s[0]==1:
3334              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3335                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3336              out[0,0]=1./arg[0,0]
3337          elif s[0]==2:
3338              A11=arg[0,0]
3339              A12=arg[0,1]
3340              A21=arg[1,0]
3341              A22=arg[1,1]
3342              D = A11*A22-A12*A21
3343              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3344                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3345              D=1./D
3346              out[0,0]= A22*D
3347              out[1,0]=-A21*D
3348              out[0,1]=-A12*D
3349              out[1,1]= A11*D
3350          elif s[0]==3:
3351              A11=arg[0,0]
3352              A21=arg[1,0]
3353              A31=arg[2,0]
3354              A12=arg[0,1]
3355              A22=arg[1,1]
3356              A32=arg[2,1]
3357              A13=arg[0,2]
3358              A23=arg[1,2]
3359              A33=arg[2,2]
3360              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3361              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3362                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3363              D=1./D
3364              out[0,0]=(A22*A33-A23*A32)*D
3365              out[1,0]=(A31*A23-A21*A33)*D
3366              out[2,0]=(A21*A32-A31*A22)*D
3367              out[0,1]=(A13*A32-A12*A33)*D
3368              out[1,1]=(A11*A33-A31*A13)*D
3369              out[2,1]=(A12*A31-A11*A32)*D
3370              out[0,2]=(A12*A23-A13*A22)*D
3371              out[1,2]=(A13*A21-A11*A23)*D
3372              out[2,2]=(A11*A22-A12*A21)*D
3373          else:
3374             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3375          return out
3376    
3377    class Inverse_Symbol(DependendSymbol):
3378       """
3379       L{Symbol} representing the result of the inverse function
3380       """
3381       def __init__(self,arg):
3382          """
3383          initialization of inverse L{Symbol} with argument arg
3384          @param arg: argument of function
3385          @type arg: L{Symbol}.
3386          """
3387          if not arg.getRank()==2:
3388            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3389          s=arg.getShape()
3390          if not s[0] == s[1]:
3391            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3392          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3393    
3394       def getMyCode(self,argstrs,format="escript"):
3395          """
3396          returns a program code that can be used to evaluate the symbol.
3397    
3398          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3399          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3400          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3401          @type format: C{str}
3402          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3403          @rtype: C{str}
3404          @raise: NotImplementedError: if the requested format is not available
3405          """
3406          if format=="escript" or format=="str"  or format=="text":
3407             return "inverse(%s)"%argstrs[0]
3408          else:
3409             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3410    
3411       def substitute(self,argvals):
3412          """
3413          assigns new values to symbols in the definition of the symbol.
3414          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3415    
3416          @param argvals: new values assigned to symbols
3417          @type argvals: C{dict} with keywords of type L{Symbol}.
3418          @return: result of the substitution process. Operations are executed as much as possible.
3419          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3420          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3421          """
3422          if argvals.has_key(self):
3423             arg=argvals[self]
3424             if self.isAppropriateValue(arg):
3425                return arg
3426             else:
3427                raise TypeError,"%s: new value is not appropriate."%str(self)
3428          else:
3429             arg=self.getSubstitutedArguments(argvals)
3430             return inverse(arg[0])
3431    
3432       def diff(self,arg):
3433          """
3434          differential of this object
3435    
3436          @param arg: the derivative is calculated with respect to arg
3437          @type arg: L{escript.Symbol}
3438          @return: derivative with respect to C{arg}
3439          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3440          """
3441          if arg==self:
3442             return identity(self.getShape())
3443          else:
3444             return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)
3445    
3446    def eigenvalues(arg):
3447        """
3448        returns the eigenvalues of the square matrix arg.
3449    
3450        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3451                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3452        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3453        @return: the eigenvalues in increasing order.
3454        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3455        @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3456        """
3457        if isinstance(arg,numarray.NumArray):
3458          if not arg.rank==2:
3459            raise ValueError,"eigenvalues: argument must have rank 2"
3460          s=arg.shape      
3461          if not s[0] == s[1]:
3462            raise ValueError,"eigenvalues: argument must be a square matrix."
3463          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3464          out.sort()
3465          return out
3466        elif isinstance(arg,escript.Data):
3467          return escript_eigenvalues(arg)
3468        elif isinstance(arg,Symbol):
3469          if not arg.getRank()==2:
3470            raise ValueError,"eigenvalues: argument must have rank 2"
3471          s=arg.getShape()      
3472          if not s[0] == s[1]:
3473            raise ValueError,"eigenvalues: argument must be a square matrix."
3474          if s[0]==1:
3475              return arg[0]
3476          elif s[0]==2:
3477              A11=arg[0,0]
3478              A12=arg[0,1]
3479              A22=arg[1,1]
3480              trA=(A11+A22)/2.
3481              A11-=trA
3482              A22-=trA
3483              s=sqrt(A12**2-A11*A22)
3484              return trA+s*numarray.array([-1.,1.])
3485          elif s[0]==3:
3486              A11=arg[0,0]
3487              A12=arg[0,1]
3488              A22=arg[1,1]
3489              A13=arg[0,2]
3490              A23=arg[1,2]
3491              A33=arg[2,2]
3492              trA=(A11+A22+A33)/3.
3493              A11-=trA
3494              A22-=trA
3495              A33-=trA
3496              A13_2=A13**2
3497              A23_2=A23**2
3498              A12_2=A12**2
3499              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3500              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3501              sq_p=sqrt(p/3.)
3502              alpha_3=acos(-q*sq_p**(-3.)/2.)/3.
3503              sq_p*=2.
3504              f=cos(alpha_3)               *numarray.array([0.,0.,1.]) \
3505               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.]) \
3506               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.])
3507              return trA+sq_p*f
3508          else:
3509             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3510        elif isinstance(arg,float):
3511          return arg
3512        elif isinstance(arg,int):
3513          return float(arg)
3514        else:
3515          raise TypeError,"eigenvalues: Unknown argument type."
3516    
3517    def escript_eigenvalues(arg): # this should be implemented in C++ arg and LAPACK is data object
3518          if not arg.getRank()==2:
3519            raise ValueError,"eigenvalues: argument must have rank 2"
3520          s=arg.getShape()      
3521          if not s[0] == s[1]:
3522            raise ValueError,"eigenvalues: argument must be a square matrix."
3523          if s[0]==1:
3524              return arg[0]
3525          elif s[0]==2:
3526              A11=arg[0,0]
3527              A12=arg[0,1]
3528              A22=arg[1,1]
3529              trA=(A11+A22)/2.
3530              A11-=trA
3531              A22-=trA
3532              s=sqrt(A12**2-A11*A22)
3533              return trA+s*numarray.array([-1.,1.])
3534          elif s[0]==3:
3535              A11=arg[0,0]
3536              A12=arg[0,1]
3537              A22=arg[1,1]
3538              A13=arg[0,2]
3539              A23=arg[1,2]
3540              A33=arg[2,2]
3541              trA=(A11+A22+A33)/3.
3542              A11-=trA
3543              A22-=trA
3544              A33-=trA
3545              A13_2=A13**2
3546              A23_2=A23**2
3547              A12_2=A12**2
3548              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3549              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3550              sq_p=sqrt(p/3.)
3551              alpha_3=acos(-q*sq_p**(-3.)/2.)/3.
3552              sq_p*=2.
3553              f=escript.Data(0.,(3,),arg.getFunctionSpace())
3554              f[0]=-cos(alpha_3-numarray.pi/3.)
3555              f[1]=-cos(alpha_3+numarray.pi/3.)
3556              f[2]=cos(alpha_3)
3557              return trA+sq_p*f
3558          else:
3559             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3560  #=======================================================  #=======================================================
3561  #  Binary operations:  #  Binary operations:
3562  #=======================================================  #=======================================================
# Line 3319  def clip(arg,minval=0.,maxval=1.): Line 4014  def clip(arg,minval=0.,maxval=1.):
4014      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4015      @param minval: lower range      @param minval: lower range
4016      @type arg: C{float}      @type arg: C{float}
4017      @param maxval: uper range      @param maxval: upper range
4018      @type arg: C{float}      @type arg: C{float}
4019      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4020               less then maxval are unchanged.               less then maxval are unchanged.
4021      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4022        @raise ValueError: if minval>maxval
4023      """      """
4024      if minval>maxval:      if minval>maxval:
4025         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
# Line 3352  def inner(arg0,arg1): Line 4048  def inner(arg0,arg1):
4048      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4049      if not sh0==sh1:      if not sh0==sh1:
4050          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4051      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4052    
4053  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4054      """      """
# Line 3380  def matrixmult(arg0,arg1): Line 4076  def matrixmult(arg0,arg1):
4076          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4077      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4078          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4079      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4080    
4081  def outer(arg0,arg1):  def outer(arg0,arg1):
4082      """      """
# Line 3398  def outer(arg0,arg1): Line 4094  def outer(arg0,arg1):
4094      @return: the outer product of arg0 and arg1 at each data point      @return: the outer product of arg0 and arg1 at each data point
4095      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4096      """      """
4097      return generalTensorProduct(arg0,arg1,offset=0)      return generalTensorProduct(arg0,arg1,axis_offset=0)
4098    
4099    
4100  def tensormult(arg0,arg1):  def tensormult(arg0,arg1):
# Line 3440  def tensormult(arg0,arg1): Line 4136  def tensormult(arg0,arg1):
4136      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
4137      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4138      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4139         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4140      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4141         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4142      else:      else:
4143          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensormult: first argument must have rank 2 or 4"
4144    
4145  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4146      """      """
4147      generalized tensor product      generalized tensor product
4148    
4149      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4150    
4151      where s runs through arg0.Shape[:arg0.Rank-offset]      where s runs through arg0.Shape[:arg0.Rank-axis_offset]
4152            r runs trough arg0.Shape[:offset]            r runs trough arg0.Shape[:axis_offset]
4153            t runs through arg1.Shape[offset:]            t runs through arg1.Shape[axis_offset:]
4154    
4155      In the first case the the second dimension of arg0 and the length of arg1 must match and        In the first case the the second dimension of arg0 and the length of arg1 must match and  
4156      in the second case the two last dimensions of arg0 must match the shape of arg1.      in the second case the two last dimensions of arg0 must match the shape of arg1.
# Line 3471  def generalTensorProduct(arg0,arg1,offse Line 4167  def generalTensorProduct(arg0,arg1,offse
4167      # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols      # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4168      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4169         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4170             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4171         else:         else:
4172             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4173                 raise ValueError,"generalTensorProduct: dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)                 raise ValueError,"generalTensorProduct: dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4174             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4175             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4176             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4177             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4178             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4179             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4180             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4181             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4182             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4183             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float)
4184             for i0 in range(d0):             for i0 in range(d0):
4185                      for i1 in range(d1):                      for i1 in range(d1):
4186                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4187             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4188             return out             return out
4189      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4190         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4191             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4192         else:         else:
4193             return escript_generalTensorProduct(arg0,arg1,offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,offset)             return escript_generalTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4194      else:            else:      
4195         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4196                                    
4197  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4198     """     """
4199     Symbol representing the quotient of two arguments.     Symbol representing the quotient of two arguments.
4200     """     """
4201     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4202         """         """
4203         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the quotient of two arguments
4204    
# Line 3515  class GeneralTensorProduct_Symbol(Depend Line 4211  class GeneralTensorProduct_Symbol(Depend
4211         """         """
4212         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
4213         sh_arg1=pokeShape(arg1)         sh_arg1=pokeShape(arg1)
4214         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4215         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4216         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4217         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4218         if not sh01==sh10:         if not sh01==sh10:
4219             raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)             raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4220         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,offset])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4221    
4222     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4223        """        """
# Line 3536  class GeneralTensorProduct_Symbol(Depend Line 4232  class GeneralTensorProduct_Symbol(Depend
4232        @raise: NotImplementedError: if the requested format is not available        @raise: NotImplementedError: if the requested format is not available
4233        """        """
4234        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4235           return "generalTensorProduct(%s,%s,offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4236        else:        else:
4237           raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)           raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4238    
# Line 3561  class GeneralTensorProduct_Symbol(Depend Line 4257  class GeneralTensorProduct_Symbol(Depend
4257           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4258           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4259    
4260  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct
4261      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4262      # calculate the return shape:      # calculate the return shape:
4263      shape0=arg0.getShape()[:arg0.getRank()-offset]      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4264      shape01=arg0.getShape()[arg0.getRank()-offset:]      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4265      shape10=arg1.getShape()[:offset]      shape10=arg1.getShape()[:axis_offset]
4266      shape1=arg1.getShape()[offset:]      shape1=arg1.getShape()[axis_offset:]
4267      if not shape01==shape10:      if not shape01==shape10:
4268          raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)          raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4269    
4270        # whatr function space should be used? (this here is not good!)
4271        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4272      # create return value:      # create return value:
4273      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())      out=escript.Data(0.,tuple(shape0+shape1),fs)
4274      #      #
4275      s0=[[]]      s0=[[]]
4276      for k in shape0:      for k in shape0:
# Line 3595  def escript_generalTensorProduct(arg0,ar Line 4293  def escript_generalTensorProduct(arg0,ar
4293    
4294      for i0 in s0:      for i0 in s0:
4295         for i1 in s1:         for i1 in s1:
4296           s=escript.Scalar(0.,arg0.getFunctionSpace())           s=escript.Scalar(0.,fs)
4297           for i01 in s01:           for i01 in s01:
4298              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))
4299           out.__setitem__(tuple(i0+i1),s)           out.__setitem__(tuple(i0+i1),s)
4300      return out      return out
4301    
4302    
4303  #=========================================================  #=========================================================
4304  #   some little helpers  #  functions dealing with spatial dependency
4305  #=========================================================  #=========================================================
4306  def grad(arg,where=None):  def grad(arg,where=None):
4307      """      """
4308      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4309    
4310      @param arg:   Data object representing the function which gradient      If C{g} is the returned object, then
4311                    to be calculated.  
4312          - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4313          - if C{arg} is rank 1 C{g[i,s]} is the derivative of C{arg[i]} with respect to the C{s}-th spatial dimension.
4314          - if C{arg} is rank 2 C{g[i,j,s]} is the derivative of C{arg[i,j]} with respect to the C{s}-th spatial dimension.
4315          - if C{arg} is rank 3 C{g[i,j,k,s]} is the derivative of C{arg[i,j,k]} with respect to the C{s}-th spatial dimension.
4316    
4317        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4318        @type arg: L{escript.Data} or L{Symbol}
4319      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4320                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4321        @type where: C{None} or L{escript.FunctionSpace}
4322        @return: gradient of arg.
4323        @rtype:  L{escript.Data} or L{Symbol}
4324      """      """
4325      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4326         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3621  def grad(arg,where=None): Line 4330  def grad(arg,where=None):
4330         else:         else:
4331            return arg._grad(where)            return arg._grad(where)
4332      else:      else:
4333        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4334    
4335    class Grad_Symbol(DependendSymbol):
4336       """
4337       L{Symbol} representing the result of the gradient operator
4338       """
4339       def __init__(self,arg,where=None):
4340          """
4341          initialization of gradient L{Symbol} with argument arg
4342          @param arg: argument of function
4343          @type arg: L{Symbol}.
4344          @param where: FunctionSpace in which the gradient will be calculated.
4345                      If not present or C{None} an appropriate default is used.
4346          @type where: C{None} or L{escript.FunctionSpace}
4347          """
4348          d=arg.getDim()
4349          if d==None:
4350             raise ValueError,"argument must have a spatial dimension"
4351          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4352    
4353       def getMyCode(self,argstrs,format="escript"):
4354          """
4355          returns a program code that can be used to evaluate the symbol.
4356    
4357          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4358          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4359          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4360          @type format: C{str}
4361          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4362          @rtype: C{str}
4363          @raise: NotImplementedError: if the requested format is not available
4364          """
4365          if format=="escript" or format=="str"  or format=="text":
4366             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4367          else:
4368             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4369    
4370       def substitute(self,argvals):
4371          """
4372          assigns new values to symbols in the definition of the symbol.
4373          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4374    
4375          @param argvals: new values assigned to symbols
4376          @type argvals: C{dict} with keywords of type L{Symbol}.
4377          @return: result of the substitution process. Operations are executed as much as possible.
4378          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4379          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4380          """
4381          if argvals.has_key(self):
4382             arg=argvals[self]
4383             if self.isAppropriateValue(arg):
4384                return arg
4385             else:
4386                raise TypeError,"%s: new value is not appropriate."%str(self)
4387          else:
4388             arg=self.getSubstitutedArguments(argvals)
4389             return grad(arg[0],where=arg[1])
4390    
4391       def diff(self,arg):
4392          """
4393          differential of this object
4394    
4395          @param arg: the derivative is calculated with respect to arg
4396          @type arg: L{escript.Symbol}
4397          @return: derivative with respect to C{arg}
4398          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4399          """
4400          if arg==self:
4401             return identity(self.getShape())
4402          else:
4403             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4404    
4405  def integrate(arg,where=None):  def integrate(arg,where=None):
4406      """      """
4407      Return the integral if the function represented by Data object arg over      Return the integral of the function C{arg} over its domain. If C{where} is present C{arg} is interpolated to C{where}
4408      its domain.      before integration.
4409    
4410      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4411        @type arg: L{escript.Data} or L{Symbol}
4412      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4413                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4414        @type where: C{None} or L{escript.FunctionSpace}
4415        @return: integral of arg.
4416        @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}
4417      """      """
4418      if isinstance(arg,numarray.NumArray):      if isinstance(arg,Symbol):
         if checkForZero(arg):  
            return arg  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,float):  
         if checkForZero(arg):  
            return arg  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,int):  
         if checkForZero(arg):  
            return float(arg)  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,Symbol):  
4419         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4420      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4421         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
# Line 3658  def integrate(arg,where=None): Line 4426  def integrate(arg,where=None):
4426      else:      else:
4427        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4428    
4429    class Integrate_Symbol(DependendSymbol):
4430       """
4431       L{Symbol} representing the result of the spatial integration operator
4432       """
4433       def __init__(self,arg,where=None):
4434          """
4435          initialization of integration L{Symbol} with argument arg
4436          @param arg: argument of the integration
4437          @type arg: L{Symbol}.
4438          @param where: FunctionSpace in which the integration will be calculated.
4439                      If not present or C{None} an appropriate default is used.
4440          @type where: C{None} or L{escript.FunctionSpace}
4441          """
4442          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4443    
4444       def getMyCode(self,argstrs,format="escript"):
4445          """
4446          returns a program code that can be used to evaluate the symbol.
4447    
4448          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4449          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4450          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4451          @type format: C{str}
4452          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4453          @rtype: C{str}
4454          @raise: NotImplementedError: if the requested format is not available
4455          """
4456          if format=="escript" or format=="str"  or format=="text":
4457             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4458          else:
4459             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4460    
4461       def substitute(self,argvals):
4462          """
4463          assigns new values to symbols in the definition of the symbol.
4464          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4465    
4466          @param argvals: new values assigned to symbols
4467          @type argvals: C{dict} with keywords of type L{Symbol}.
4468          @return: result of the substitution process. Operations are executed as much as possible.
4469          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4470          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4471          """
4472          if argvals.has_key(self):
4473             arg=argvals[self]
4474             if self.isAppropriateValue(arg):
4475                return arg
4476             else:
4477                raise TypeError,"%s: new value is not appropriate."%str(self)
4478          else:
4479             arg=self.getSubstitutedArguments(argvals)
4480             return integrate(arg[0],where=arg[1])
4481    
4482       def diff(self,arg):
4483          """
4484          differential of this object
4485    
4486          @param arg: the derivative is calculated with respect to arg
4487          @type arg: L{escript.Symbol}
4488          @return: derivative with respect to C{arg}
4489          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4490          """
4491          if arg==self:
4492             return identity(self.getShape())
4493          else:
4494             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4495    
4496    
4497  def interpolate(arg,where):  def interpolate(arg,where):
4498      """      """
4499      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
4500    
4501      @param arg:    interpolant      @param arg: interpolant
4502      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
4503        @param where: FunctionSpace to be interpolated to
4504        @type where: L{escript.FunctionSpace}
4505        @return: interpolated argument
4506        @rtype:  C{escript.Data} or L{Symbol}
4507      """      """
4508      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4509         return Interpolated_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
4510      else:      else:
4511         return escript.Data(arg,where)         return escript.Data(arg,where)
4512    
4513  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
4514      """     """
4515      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
4516       """
4517       def __init__(self,arg,where):
4518          """
4519          initialization of interpolation L{Symbol} with argument arg
4520          @param arg: argument of the interpolation
4521          @type arg: L{Symbol}.
4522          @param where: FunctionSpace into which the argument is interpolated.
4523          @type where: L{escript.FunctionSpace}
4524          """
4525          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4526    
4527      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
4528                    be calculated.        """
4529      @param where: FunctionSpace in which the gradient will be calculated.        returns a program code that can be used to evaluate the symbol.
                   If not present or C{None} an appropriate default is used.  
     """  
     g=grad(arg,where)  
     return trace(g,axis0=g.getRank()-2,axis1=g.getRank()-1)  
4530    
4531  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
4532      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4533      Returns the jump of arg across a continuity.        @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4534          @type format: C{str}
4535          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4536          @rtype: C{str}
4537          @raise: NotImplementedError: if the requested format is not available
4538          """
4539          if format=="escript" or format=="str"  or format=="text":
4540             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
4541          else:
4542             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4543    
4544      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
4545                    to be calculated.        """
4546      """        assigns new values to symbols in the definition of the symbol.
4547      d=arg.getDomain()        The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
     return arg.interpolate(escript.FunctionOnContactOne(d))-arg.interpolate(escript.FunctionOnContactZero(d))  
4548    
4549  #=============================        @param argvals: new values assigned to symbols
4550  #        @type argvals: C{dict} with keywords of type L{Symbol}.
4551  # wrapper for various functions: if the argument has attribute the function name        @return: result of the substitution process. Operations are executed as much as possible.
4552  # as an argument it calls the corresponding methods. Otherwise the corresponding        @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4553  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4554          """
4555          if argvals.has_key(self):
4556             arg=argvals[self]
4557             if self.isAppropriateValue(arg):
4558                return arg
4559             else:
4560                raise TypeError,"%s: new value is not appropriate."%str(self)
4561          else:
4562             arg=self.getSubstitutedArguments(argvals)
4563             return interpolate(arg[0],where=arg[1])
4564    
4565  # functions involving the underlying Domain:     def diff(self,arg):
4566          """
4567          differential of this object
4568    
4569  def transpose(arg,axis=None):        @param arg: the derivative is calculated with respect to arg
4570          @type arg: L{escript.Symbol}
4571          @return: derivative with respect to C{arg}
4572          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
4573          """
4574          if arg==self:
4575             return identity(self.getShape())
4576          else:
4577             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4578    
4579    
4580    def div(arg,where=None):
4581      """      """
4582      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
4583    
4584      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
4585        @type arg: L{escript.Data} or L{Symbol}
4586        @param where: FunctionSpace in which the divergence will be calculated.
4587                      If not present or C{None} an appropriate default is used.
4588        @type where: C{None} or L{escript.FunctionSpace}
4589        @return: divergence of arg.
4590        @rtype:  L{escript.Data} or L{Symbol}
4591      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
4592      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4593         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
4594      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4595         # hack for transpose          dim=arg.getDomain().getDim()
        r=arg.getRank()  
        if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"  
        s=arg.getShape()  
        out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())  
        for i in range(s[0]):  
           for j in range(s[1]):  
              out[j,i]=arg[i,j]  
        return out  
        # end hack for transpose  
        return arg.transpose(axis)  
4596      else:      else:
4597         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
4598        if not arg.getShape()==(dim,):
4599          raise ValueError,"div: expected shape is (%s,)"%dim
4600        return trace(grad(arg,where))
4601    
4602  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
4603      """      """
4604      Return      returns the jump of arg across the continuity of the domain
4605    
4606      @param arg:      @param arg: argument
4607        @type arg: L{escript.Data} or L{Symbol}
4608        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
4609                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
4610        @type domain: C{None} or L{escript.Domain}
4611        @return: jump of arg
4612        @rtype:  L{escript.Data} or L{Symbol}
4613      """      """
4614      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
4615         s=list(arg.getShape())              return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
        s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])  
        return Trace_Symbol(arg,axis0=axis0,axis1=axis1)  
     elif isinstance(arg,escript.Data):  
        # hack for trace  
        s=arg.getShape()  
        if s[axis0]!=s[axis1]:  
            raise ValueError,"illegal axis in trace"  
        out=escript.Scalar(0.,arg.getFunctionSpace())  
        for i in range(s[axis0]):  
           out+=arg[i,i]  
        return out  
        # end hack for trace  
     else:  
        return numarray.trace(arg,axis0=axis0,axis1=axis1)  
4616    
4617    def L2(arg):
4618        """
4619        returns the L2 norm of arg at where
4620        
4621        @param arg: function which L2 to be calculated.
4622        @type arg: L{escript.Data} or L{Symbol}
4623        @return: L2 norm of arg.
4624        @rtype:  L{float} or L{Symbol}
4625        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
4626        """
4627        return sqrt(integrate(inner(arg,arg)))
4628    #=============================
4629    #
4630    
4631  def reorderComponents(arg,index):  def reorderComponents(arg,index):
4632      """      """
4633      resorts the component of arg according to index      resorts the component of arg according to index
4634    
4635      """      """
4636      pass      raise NotImplementedError
4637  #  #
4638  # $Log: util.py,v $  # $Log: util.py,v $
4639  # Revision 1.14.2.16  2005/10/19 06:09:57  gross  # Revision 1.14.2.16  2005/10/19 06:09:57  gross

Legend:
Removed from v.396  
changed lines
  Added in v.536

  ViewVC Help
Powered by ViewVC 1.1.26