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

revision 429 by gross, Wed Jan 11 05:53:40 2006 UTC revision 528 by gross, Wed Feb 15 02:20:33 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 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):

49  #  #
50  # slicing: get  # slicing: get
51  #          set  #          set
# Line 122  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 144  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 160  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 175  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 188  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 830  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 880  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 3013  class Trace_Symbol(DependendSymbol): Line 3110  class Trace_Symbol(DependendSymbol):
3110        else:        else:
3111           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           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
3281    def inverse(arg):
3282        """
3283        returns the inverse of the square matrix arg.
3284
3285        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3286        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3287        @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3288        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3289        @remark: for L{escript.Data} objects the dimension is restricted to 3.
3290        """
3291        if isinstance(arg,numarray.NumArray):
3292          return numarray.linear_algebra.inverse(arg)
3293        elif isinstance(arg,escript.Data):
3294          return escript_inverse(arg)
3295        elif isinstance(arg,float):
3296          return 1./arg
3297        elif isinstance(arg,int):
3298          return 1./float(arg)
3299        elif isinstance(arg,Symbol):
3300          return Inverse_Symbol(arg)
3301        else:
3302          raise TypeError,"inverse: Unknown argument type."
3303
3304    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3305          "arg is a Data objects!!!"
3306          if not arg.getRank()==2:
3307            raise ValueError,"escript_inverse: argument must have rank 2"
3308          s=arg.getShape()
3309          if not s[0] == s[1]:
3310            raise ValueError,"escript_inverse: argument must be a square matrix."
3311          out=escript.Data(0.,s,arg.getFunctionSpace())
3312          if s[0]==1:
3313              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3314                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3315              out[0,0]=1./arg[0,0]
3316          elif s[0]==2:
3317              A11=arg[0,0]
3318              A12=arg[0,1]
3319              A21=arg[1,0]
3320              A22=arg[1,1]
3321              D = A11*A22-A12*A21
3322              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3323                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3324              D=1./D
3325              out[0,0]= A22*D
3326              out[1,0]=-A21*D
3327              out[0,1]=-A12*D
3328              out[1,1]= A11*D
3329          elif s[0]==3:
3330              A11=arg[0,0]
3331              A21=arg[1,0]
3332              A31=arg[2,0]
3333              A12=arg[0,1]
3334              A22=arg[1,1]
3335              A32=arg[2,1]
3336              A13=arg[0,2]
3337              A23=arg[1,2]
3338              A33=arg[2,2]
3339              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3340              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3341                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3342              D=1./D
3343              out[0,0]=(A22*A33-A23*A32)*D
3344              out[1,0]=(A31*A23-A21*A33)*D
3345              out[2,0]=(A21*A32-A31*A22)*D
3346              out[0,1]=(A13*A32-A12*A33)*D
3347              out[1,1]=(A11*A33-A31*A13)*D
3348              out[2,1]=(A12*A31-A11*A32)*D
3349              out[0,2]=(A12*A23-A13*A22)*D
3350              out[1,2]=(A13*A21-A11*A23)*D
3351              out[2,2]=(A11*A22-A12*A21)*D
3352          else:
3353             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3354          return out
3355
3356    class Inverse_Symbol(DependendSymbol):
3357       """
3358       L{Symbol} representing the result of the inverse function
3359       """
3360       def __init__(self,arg):
3361          """
3362          initialization of inverse L{Symbol} with argument arg
3363          @param arg: argument of function
3364          @type arg: L{Symbol}.
3365          """
3366          if not arg.getRank()==2:
3367            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3368          s=arg.getShape()
3369          if not s[0] == s[1]:
3370            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3371          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3372
3373       def getMyCode(self,argstrs,format="escript"):
3374          """
3375          returns a program code that can be used to evaluate the symbol.
3376
3377          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3378          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3379          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3380          @type format: C{str}
3381          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3382          @rtype: C{str}
3383          @raise: NotImplementedError: if the requested format is not available
3384          """
3385          if format=="escript" or format=="str"  or format=="text":
3386             return "inverse(%s)"%argstrs[0]
3387          else:
3388             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3389
3390       def substitute(self,argvals):
3391          """
3392          assigns new values to symbols in the definition of the symbol.
3393          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3394
3395          @param argvals: new values assigned to symbols
3396          @type argvals: C{dict} with keywords of type L{Symbol}.
3397          @return: result of the substitution process. Operations are executed as much as possible.
3398          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3399          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3400          """
3401          if argvals.has_key(self):
3402             arg=argvals[self]
3403             if self.isAppropriateValue(arg):
3404                return arg
3405             else:
3406                raise TypeError,"%s: new value is not appropriate."%str(self)
3407          else:
3408             arg=self.getSubstitutedArguments(argvals)
3409             return inverse(arg[0])
3410
3411       def diff(self,arg):
3412          """
3413          differential of this object
3414
3415          @param arg: the derivative is calculated with respect to arg
3416          @type arg: L{escript.Symbol}
3417          @return: derivative with respect to C{arg}
3418          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3419          """
3420          if arg==self:
3421             return identity(self.getShape())
3422          else:
3423             return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)
3424
3425    def eigenvalues(arg):
3426        """
3427        returns the eigenvalues of the square matrix arg.
3428
3429        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3430                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3431        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol},L {float} or L{int}
3432        @return: the list of the eigenvalues in increasing order.
3433        @rtype: C{list} of L{float}, L{escript.Data}, L{Symbol} depending on the input.
3434        @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3435        """
3436        if isinstance(arg,numarray.NumArray):
3437          if not arg.rank==2:
3438            raise ValueError,"eigenvalues: argument must have rank 2"
3439          s=arg.shape
3440          if not s[0] == s[1]:
3441            raise ValueError,"eigenvalues: argument must be a square matrix."
3442          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3443          out.sort()
3444          return numarray.array2list(out)
3445        elif isinstance(arg,escript.Data) or isinstance(arg,Symbol):
3446          if not arg.getRank()==2:
3447            raise ValueError,"eigenvalues: argument must have rank 2"
3448          s=arg.getShape()
3449          if not s[0] == s[1]:
3450            raise ValueError,"eigenvalues: argument must be a square matrix."
3451          if s[0]==1:
3452              return [arg[0,0]]
3453          elif s[0]==2:
3454              A11=arg[0,0]
3455              A12=arg[0,1]
3456              A22=arg[1,1]
3457              trA=(A11+A22)/2.
3458              A11-=trA
3459              A22-=trA
3460              s=sqrt(A12**2-A11*A22)
3461              return [trA-s,trA+s]
3462          elif s[0]==3:
3463              A11=arg[0,0]
3464              A12=arg[0,1]
3465              A22=arg[1,1]
3466              A13=arg[0,2]
3467              A23=arg[1,2]
3468              A33=arg[2,2]
3469              trA=(A11+A22+A22)/3.
3470              A11-=trA
3471              A22-=trA
3472              A33-=trA
3473              A13_2=A13**2
3474              A23_2=A23**2
3475              A12_2=A12**2
3476              p=A13_2+A_23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3477              q=A13_2*A22+A_23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3478              sq_p=sqrt(p/3.)
3479              alpha_3=acos(-q/sq_p**(1./3.)/2.)/3.
3480              sq_p*=2.
3481              return [trA+sq_p*cos(alpha_3),trA-sq_p*cos(alpha_3+numarray.pi/3.),trA-sq_p*cos(alpha_3-numarray.pi/3.)]
3482          else:
3483             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3484        elif isinstance(arg,float):
3485          return [arg]
3486        elif isinstance(arg,int):
3487          return [float(arg)]
3488        else:
3489          raise TypeError,"eigenvalues: Unknown argument type."
3490  #=======================================================  #=======================================================
3491  #  Binary operations:  #  Binary operations:
3492  #=======================================================  #=======================================================
4278        d=arg.getDim()        d=arg.getDim()
4279        if d==None:        if d==None:
4280           raise ValueError,"argument must have a spatial dimension"           raise ValueError,"argument must have a spatial dimension"
4282
4283     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4284        """        """
# Line 4045  def div(arg,where=None): Line 4519  def div(arg,where=None):
4519      @return: divergence of arg.      @return: divergence of arg.
4520      @rtype:  L{escript.Data} or L{Symbol}      @rtype:  L{escript.Data} or L{Symbol}
4521      """      """
4522      if not arg.getShape()==(arg.getDim(),):      if isinstance(arg,Symbol):
4523        raise ValueError,"div: expected shape is (%s,)"%arg.getDim()          dim=arg.getDim()
4524        elif isinstance(arg,escript.Data):
4525            dim=arg.getDomain().getDim()
4526        else:
4527            raise TypeError,"div: argument type not supported"
4528        if not arg.getShape()==(dim,):
4529          raise ValueError,"div: expected shape is (%s,)"%dim
4531
4532  def jump(arg,domain=None):  def jump(arg,domain=None):
# Line 4063  def jump(arg,domain=None): Line 4543  def jump(arg,domain=None):
4543      """      """
4544      if domain==None: domain=arg.getDomain()      if domain==None: domain=arg.getDomain()
4545      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
#=============================
#
# wrapper for various functions: if the argument has attribute the function name
# as an argument it calls the corresponding methods. Otherwise the corresponding
# numarray function is called.
4546
4547  # functions involving the underlying Domain:  def L2(arg):

def transpose(arg,axis=None):
4548      """      """
4549      Returns the transpose of the Data object arg.      returns the L2 norm of arg at where
4550
4551      @param arg:      @param arg: function which L2 to be calculated.
4552        @type arg: L{escript.Data} or L{Symbol}
4553        @return: L2 norm of arg.
4554        @rtype:  L{float} or L{Symbol}
4555        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
4556      """      """
4557      if axis==None:      return sqrt(integrate(inner(arg,arg)))
4558         r=0  #=============================
4559         if hasattr(arg,"getRank"): r=arg.getRank()  #
if hasattr(arg,"rank"): r=arg.rank
axis=r/2
if isinstance(arg,Symbol):
return Transpose_Symbol(arg,axis=r)
if isinstance(arg,escript.Data):
# hack for transpose
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)
else:
return numarray.transpose(arg,axis=axis)

4560
4561  def reorderComponents(arg,index):  def reorderComponents(arg,index):
4562      """      """
4563      resorts the component of arg according to index      resorts the component of arg according to index
4564
4565      """      """
4566      pass      raise NotImplementedError
4567  #  #
4568  # \$Log: util.py,v \$  # \$Log: util.py,v \$
4569  # 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.429 changed lines Added in v.528