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

revision 290 by gross, Fri Dec 2 03:09:32 2005 UTC revision 587 by gross, Fri Mar 10 02:26:50 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 minval(arg0):
# def maxval(arg0):
# def transpose(arg,axis=None):
# def trace(arg,axis0=0,axis1=1):
# def length(arg):
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 127  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 145  def identity(shape=()): Line 133  def identity(shape=()):
133     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
134     """     """
135     if len(shape)>0:     if len(shape)>0:
136        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
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 165  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 180  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 193  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 365  def testForZero(arg): Line 358  def testForZero(arg):
358      @return : True if the argument is identical to zero.      @return : True if the argument is identical to zero.
359      @rtype : C{bool}      @rtype : C{bool}
360      """      """
361      try:      if isinstance(arg,numarray.NumArray):
362         return not Lsup(arg)>0.         return not Lsup(arg)>0.
363      except TypeError:      elif isinstance(arg,escript.Data):
364           return False
365        elif isinstance(arg,float):
366           return not Lsup(arg)>0.
367        elif isinstance(arg,int):
368           return not Lsup(arg)>0.
369        elif isinstance(arg,Symbol):
370           return False
371        else:
372         return False         return False
373
374  def matchType(arg0=0.,arg1=0.):  def matchType(arg0=0.,arg1=0.):
# Line 388  def matchType(arg0=0.,arg1=0.): Line 389  def matchType(arg0=0.,arg1=0.):
389         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
390            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
391         elif isinstance(arg1,float):         elif isinstance(arg1,float):
392            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
393         elif isinstance(arg1,int):         elif isinstance(arg1,int):
394            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
395         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
396            pass            pass
397         else:         else:
# Line 414  def matchType(arg0=0.,arg1=0.): Line 415  def matchType(arg0=0.,arg1=0.):
415         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
416            pass            pass
417         elif isinstance(arg1,float):         elif isinstance(arg1,float):
418            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
419         elif isinstance(arg1,int):         elif isinstance(arg1,int):
420            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
421         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
422            pass            pass
423         else:         else:
424            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."
425      elif isinstance(arg0,float):      elif isinstance(arg0,float):
426         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
427            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
428         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
429            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
430         elif isinstance(arg1,float):         elif isinstance(arg1,float):
431            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
432            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
433         elif isinstance(arg1,int):         elif isinstance(arg1,int):
434            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
435            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
436         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
437            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
438         else:         else:
439            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."
440      elif isinstance(arg0,int):      elif isinstance(arg0,int):
441         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
442            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
443         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
444            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
445         elif isinstance(arg1,float):         elif isinstance(arg1,float):
446            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
447            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
448         elif isinstance(arg1,int):         elif isinstance(arg1,int):
449            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
450            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
451         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
452            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
453         else:         else:
454            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."
455      else:      else:
# Line 471  def matchShape(arg0,arg1): Line 472  def matchShape(arg0,arg1):
472      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
473      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
474      if len(sh0)<len(sh):      if len(sh0)<len(sh):
475         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
476      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
477         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
478      else:      else:
479         return arg0,arg1         return arg0,arg1
480  #=========================================================  #=========================================================
# Line 599  class Symbol(object): Line 600  class Symbol(object):
600            else:            else:
601                s=pokeShape(s)+arg.getShape()                s=pokeShape(s)+arg.getShape()
602                if len(s)>0:                if len(s)>0:
603                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
604                else:                else:
605                   out.append(a)                   out.append(a)
606         return out         return out
# Line 689  class Symbol(object): Line 690  class Symbol(object):
690         else:         else:
691            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
692            if len(s)>0:            if len(s)>0:
693               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
694            else:            else:
695               return 0.               return 0.
696
# Line 827  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 877  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 909  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        return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1013          if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1014          return out
1015     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1016        return arg._wherePositive()        return arg._wherePositive()
1017     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 989  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        return numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1095          if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1096          return out
1097     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1098        return arg._whereNegative()        return arg._whereNegative()
1099     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1069  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        return numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1177          if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1178          return out
1179     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1180        return arg._whereNonNegative()        return arg._whereNonNegative()
1181     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1097  def whereNonPositive(arg): Line 1203  def whereNonPositive(arg):
1203     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1204     """     """
1205     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1206        return numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1207          if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1208          return out
1209     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1210        return arg._whereNonPositive()        return arg._whereNonPositive()
1211     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1127  def whereZero(arg,tol=0.): Line 1235  def whereZero(arg,tol=0.):
1235     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1236     """     """
1237     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1238        return numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1239          if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1240          return out
1241     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1242        if tol>0.:        if tol>0.:
1243           return whereNegative(abs(arg)-tol)           return whereNegative(abs(arg)-tol)
# Line 1208  def whereNonZero(arg,tol=0.): Line 1318  def whereNonZero(arg,tol=0.):
1318     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1319     """     """
1320     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1321        return numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1322          if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1323          return out
1324     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1325        if tol>0.:        if tol>0.:
1326           return 1.-whereZero(arg,tol)           return 1.-whereZero(arg,tol)
# Line 2589  def sign(arg): Line 2701  def sign(arg):
2701     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2702     """     """
2703     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
2704        return numarray.sign(arg)        return wherePositive(arg)-whereNegative(arg)
2705     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
2706        return arg._sign()        return arg._sign()
2707     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 2679  class Abs_Symbol(DependendSymbol): Line 2791  class Abs_Symbol(DependendSymbol):
2791           val=matchShape(sign(myarg),self.getDifferentiatedArguments(arg)[0])           val=matchShape(sign(myarg),self.getDifferentiatedArguments(arg)[0])
2792           return val[0]*val[1]           return val[0]*val[1]
2793
2794    def minval(arg):
2795       """
2796       returns minimum value over all components of arg at each data point
2797
2798       @param arg: argument
2799       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2800       @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2801       @raises TypeError: if the type of the argument is not expected.
2802       """
2803       if isinstance(arg,numarray.NumArray):
2804          if arg.rank==0:
2805             return float(arg)
2806          else:
2807             return arg.min()
2808       elif isinstance(arg,escript.Data):
2809          return arg._minval()
2810       elif isinstance(arg,float):
2811          return arg
2812       elif isinstance(arg,int):
2813          return float(arg)
2814       elif isinstance(arg,Symbol):
2815          return Minval_Symbol(arg)
2816       else:
2817          raise TypeError,"minval: Unknown argument type."
2818
2819    class Minval_Symbol(DependendSymbol):
2820       """
2821       L{Symbol} representing the result of the minimum value function
2822       """
2823       def __init__(self,arg):
2824          """
2825          initialization of minimum value L{Symbol} with argument arg
2826          @param arg: argument of function
2827          @type arg: typically L{Symbol}.
2828          """
2829          DependendSymbol.__init__(self,args=[arg],shape=(),dim=arg.getDim())
2830
2831       def getMyCode(self,argstrs,format="escript"):
2832          """
2833          returns a program code that can be used to evaluate the symbol.
2834
2835          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2836          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2837          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2838          @type format: C{str}
2839          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2840          @rtype: C{str}
2841          @raise: NotImplementedError: if the requested format is not available
2842          """
2843          if isinstance(argstrs,list):
2844              argstrs=argstrs[0]
2845          if format=="escript" or format=="str"  or format=="text":
2846             return "minval(%s)"%argstrs
2847          else:
2848             raise NotImplementedError,"Minval_Symbol does not provide program code for format %s."%format
2849
2850       def substitute(self,argvals):
2851          """
2852          assigns new values to symbols in the definition of the symbol.
2853          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2854
2855          @param argvals: new values assigned to symbols
2856          @type argvals: C{dict} with keywords of type L{Symbol}.
2857          @return: result of the substitution process. Operations are executed as much as possible.
2858          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2859          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2860          """
2861          if argvals.has_key(self):
2862             arg=argvals[self]
2863             if self.isAppropriateValue(arg):
2864                return arg
2865             else:
2866                raise TypeError,"%s: new value is not appropriate."%str(self)
2867          else:
2868             arg=self.getSubstitutedArguments(argvals)[0]
2869             return minval(arg)
2870
2871    def maxval(arg):
2872       """
2873       returns maximum value over all components of arg at each data point
2874
2875       @param arg: argument
2876       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2877       @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2878       @raises TypeError: if the type of the argument is not expected.
2879       """
2880       if isinstance(arg,numarray.NumArray):
2881          if arg.rank==0:
2882             return float(arg)
2883          else:
2884             return arg.max()
2885       elif isinstance(arg,escript.Data):
2886          return arg._maxval()
2887       elif isinstance(arg,float):
2888          return arg
2889       elif isinstance(arg,int):
2890          return float(arg)
2891       elif isinstance(arg,Symbol):
2892          return Maxval_Symbol(arg)
2893       else:
2894          raise TypeError,"maxval: Unknown argument type."
2895
2896    class Maxval_Symbol(DependendSymbol):
2897       """
2898       L{Symbol} representing the result of the maximum value function
2899       """
2900       def __init__(self,arg):
2901          """
2902          initialization of maximum value L{Symbol} with argument arg
2903          @param arg: argument of function
2904          @type arg: typically L{Symbol}.
2905          """
2906          DependendSymbol.__init__(self,args=[arg],shape=(),dim=arg.getDim())
2907
2908       def getMyCode(self,argstrs,format="escript"):
2909          """
2910          returns a program code that can be used to evaluate the symbol.
2911
2912          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2913          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2914          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2915          @type format: C{str}
2916          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2917          @rtype: C{str}
2918          @raise: NotImplementedError: if the requested format is not available
2919          """
2920          if isinstance(argstrs,list):
2921              argstrs=argstrs[0]
2922          if format=="escript" or format=="str"  or format=="text":
2923             return "maxval(%s)"%argstrs
2924          else:
2925             raise NotImplementedError,"Maxval_Symbol does not provide program code for format %s."%format
2926
2927       def substitute(self,argvals):
2928          """
2929          assigns new values to symbols in the definition of the symbol.
2930          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2931
2932          @param argvals: new values assigned to symbols
2933          @type argvals: C{dict} with keywords of type L{Symbol}.
2934          @return: result of the substitution process. Operations are executed as much as possible.
2935          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2936          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2937          """
2938          if argvals.has_key(self):
2939             arg=argvals[self]
2940             if self.isAppropriateValue(arg):
2941                return arg
2942             else:
2943                raise TypeError,"%s: new value is not appropriate."%str(self)
2944          else:
2945             arg=self.getSubstitutedArguments(argvals)[0]
2946             return maxval(arg)
2947
2948    def length(arg):
2949       """
2950       returns length/Euclidean norm of argument arg at each data point
2951
2952       @param arg: argument
2953       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2954       @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2955       """
2956       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.Float64)
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        if isinstance(arg,numarray.NumArray):
3290          if arg.rank==2:
3291            if not (arg.shape[0]==arg.shape[1]):
3292               raise ValueError,"symmetric: argument must be square."
3293          elif arg.rank==4:
3294            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3295               raise ValueError,"symmetric: argument must be square."
3296          else:
3297            raise ValueError,"symmetric: rank 2 or 4 is required."
3298          return (arg+transpose(arg))/2
3299        elif isinstance(arg,escript.Data):
3300          return escript_symmetric(arg)
3301        elif isinstance(arg,float):
3302          return arg
3303        elif isinstance(arg,int):
3304          return float(arg)
3305        elif isinstance(arg,Symbol):
3306          if arg.getRank()==2:
3307            if not (arg.getShape()[0]==arg.getShape()[1]):
3308               raise ValueError,"symmetric: argument must be square."
3309          elif arg.getRank()==4:
3310            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3311               raise ValueError,"symmetric: argument must be square."
3312          else:
3313            raise ValueError,"symmetric: rank 2 or 4 is required."
3314          return (arg+transpose(arg))/2
3315        else:
3316          raise TypeError,"symmetric: Unknown argument type."
3317
3318    def escript_symmetric(arg): # this should be implemented in c++
3319          if arg.getRank()==2:
3320            if not (arg.getShape()[0]==arg.getShape()[1]):
3321               raise ValueError,"escript_symmetric: argument must be square."
3322            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3323            for i0 in range(arg.getShape()[0]):
3324               for i1 in range(arg.getShape()[1]):
3325                  out[i0,i1]=(arg[i0,i1]+arg[i1,i0])/2.
3326          elif arg.getRank()==4:
3327            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3328               raise ValueError,"escript_symmetric: argument must be square."
3329            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3330            for i0 in range(arg.getShape()[0]):
3331               for i1 in range(arg.getShape()[1]):
3332                  for i2 in range(arg.getShape()[2]):
3333                     for i3 in range(arg.getShape()[3]):
3334                         out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]+arg[i2,i3,i0,i1])/2.
3335          else:
3336            raise ValueError,"escript_symmetric: rank 2 or 4 is required."
3337          return out
3338
3339    def nonsymmetric(arg):
3340        """
3341        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3342
3343        @param arg: square matrix. Must have rank 2 or 4 and be square.
3344        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3345        @return: nonsymmetric part of arg
3346        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3347        """
3348        if isinstance(arg,numarray.NumArray):
3349          if arg.rank==2:
3350            if not (arg.shape[0]==arg.shape[1]):
3351               raise ValueError,"nonsymmetric: argument must be square."
3352          elif arg.rank==4:
3353            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3354               raise ValueError,"nonsymmetric: argument must be square."
3355          else:
3356            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3357          return (arg-transpose(arg))/2
3358        elif isinstance(arg,escript.Data):
3359          return escript_nonsymmetric(arg)
3360        elif isinstance(arg,float):
3361          return arg
3362        elif isinstance(arg,int):
3363          return float(arg)
3364        elif isinstance(arg,Symbol):
3365          if arg.getRank()==2:
3366            if not (arg.getShape()[0]==arg.getShape()[1]):
3367               raise ValueError,"nonsymmetric: argument must be square."
3368          elif arg.getRank()==4:
3369            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3370               raise ValueError,"nonsymmetric: argument must be square."
3371          else:
3372            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3373          return (arg-transpose(arg))/2
3374        else:
3375          raise TypeError,"nonsymmetric: Unknown argument type."
3376
3377    def escript_nonsymmetric(arg): # this should be implemented in c++
3378          if arg.getRank()==2:
3379            if not (arg.getShape()[0]==arg.getShape()[1]):
3380               raise ValueError,"escript_nonsymmetric: argument must be square."
3381            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3382            for i0 in range(arg.getShape()[0]):
3383               for i1 in range(arg.getShape()[1]):
3384                  out[i0,i1]=(arg[i0,i1]-arg[i1,i0])/2.
3385          elif arg.getRank()==4:
3386            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3387               raise ValueError,"escript_nonsymmetric: argument must be square."
3388            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3389            for i0 in range(arg.getShape()[0]):
3390               for i1 in range(arg.getShape()[1]):
3391                  for i2 in range(arg.getShape()[2]):
3392                     for i3 in range(arg.getShape()[3]):
3393                         out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]-arg[i2,i3,i0,i1])/2.
3394          else:
3395            raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."
3396          return out
3397
3398
3399    def inverse(arg):
3400        """
3401        returns the inverse of the square matrix arg.
3402
3403        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3404        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3405        @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3406        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3407        @remark: for L{escript.Data} objects the dimension is restricted to 3.
3408        """
3409        if isinstance(arg,numarray.NumArray):
3410          return numarray.linear_algebra.inverse(arg)
3411        elif isinstance(arg,escript.Data):
3412          return escript_inverse(arg)
3413        elif isinstance(arg,float):
3414          return 1./arg
3415        elif isinstance(arg,int):
3416          return 1./float(arg)
3417        elif isinstance(arg,Symbol):
3418          return Inverse_Symbol(arg)
3419        else:
3420          raise TypeError,"inverse: Unknown argument type."
3421
3422    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3423          "arg is a Data objects!!!"
3424          if not arg.getRank()==2:
3425            raise ValueError,"escript_inverse: argument must have rank 2"
3426          s=arg.getShape()
3427          if not s[0] == s[1]:
3428            raise ValueError,"escript_inverse: argument must be a square matrix."
3429          out=escript.Data(0.,s,arg.getFunctionSpace())
3430          if s[0]==1:
3431              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3432                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3433              out[0,0]=1./arg[0,0]
3434          elif s[0]==2:
3435              A11=arg[0,0]
3436              A12=arg[0,1]
3437              A21=arg[1,0]
3438              A22=arg[1,1]
3439              D = A11*A22-A12*A21
3440              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3441                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3442              D=1./D
3443              out[0,0]= A22*D
3444              out[1,0]=-A21*D
3445              out[0,1]=-A12*D
3446              out[1,1]= A11*D
3447          elif s[0]==3:
3448              A11=arg[0,0]
3449              A21=arg[1,0]
3450              A31=arg[2,0]
3451              A12=arg[0,1]
3452              A22=arg[1,1]
3453              A32=arg[2,1]
3454              A13=arg[0,2]
3455              A23=arg[1,2]
3456              A33=arg[2,2]
3457              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3458              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3459                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3460              D=1./D
3461              out[0,0]=(A22*A33-A23*A32)*D
3462              out[1,0]=(A31*A23-A21*A33)*D
3463              out[2,0]=(A21*A32-A31*A22)*D
3464              out[0,1]=(A13*A32-A12*A33)*D
3465              out[1,1]=(A11*A33-A31*A13)*D
3466              out[2,1]=(A12*A31-A11*A32)*D
3467              out[0,2]=(A12*A23-A13*A22)*D
3468              out[1,2]=(A13*A21-A11*A23)*D
3469              out[2,2]=(A11*A22-A12*A21)*D
3470          else:
3471             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3472          return out
3473
3474    class Inverse_Symbol(DependendSymbol):
3475       """
3476       L{Symbol} representing the result of the inverse function
3477       """
3478       def __init__(self,arg):
3479          """
3480          initialization of inverse L{Symbol} with argument arg
3481          @param arg: argument of function
3482          @type arg: L{Symbol}.
3483          """
3484          if not arg.getRank()==2:
3485            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3486          s=arg.getShape()
3487          if not s[0] == s[1]:
3488            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3489          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3490
3491       def getMyCode(self,argstrs,format="escript"):
3492          """
3493          returns a program code that can be used to evaluate the symbol.
3494
3495          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3496          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3497          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3498          @type format: C{str}
3499          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3500          @rtype: C{str}
3501          @raise: NotImplementedError: if the requested format is not available
3502          """
3503          if format=="escript" or format=="str"  or format=="text":
3504             return "inverse(%s)"%argstrs[0]
3505          else:
3506             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3507
3508       def substitute(self,argvals):
3509          """
3510          assigns new values to symbols in the definition of the symbol.
3511          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3512
3513          @param argvals: new values assigned to symbols
3514          @type argvals: C{dict} with keywords of type L{Symbol}.
3515          @return: result of the substitution process. Operations are executed as much as possible.
3516          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3517          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3518          """
3519          if argvals.has_key(self):
3520             arg=argvals[self]
3521             if self.isAppropriateValue(arg):
3522                return arg
3523             else:
3524                raise TypeError,"%s: new value is not appropriate."%str(self)
3525          else:
3526             arg=self.getSubstitutedArguments(argvals)
3527             return inverse(arg[0])
3528
3529       def diff(self,arg):
3530          """
3531          differential of this object
3532
3533          @param arg: the derivative is calculated with respect to arg
3534          @type arg: L{escript.Symbol}
3535          @return: derivative with respect to C{arg}
3536          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3537          """
3538          if arg==self:
3539             return identity(self.getShape())
3540          else:
3541             return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)
3542
3543    def eigenvalues(arg):
3544        """
3545        returns the eigenvalues of the square matrix arg.
3546
3547        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3548                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3549        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3550        @return: the eigenvalues in increasing order.
3551        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3552        @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3553        """
3554        if isinstance(arg,numarray.NumArray):
3555          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3556          out.sort()
3557          return out
3558        elif isinstance(arg,escript.Data):
3559          return arg._eigenvalues()
3560        elif isinstance(arg,Symbol):
3561          if not arg.getRank()==2:
3562            raise ValueError,"eigenvalues: argument must have rank 2"
3563          s=arg.getShape()
3564          if not s[0] == s[1]:
3565            raise ValueError,"eigenvalues: argument must be a square matrix."
3566          if s[0]==1:
3567              return arg[0]
3568          elif s[0]==2:
3569              arg1=symmetric(arg)
3570              A11=arg1[0,0]
3571              A12=arg1[0,1]
3572              A22=arg1[1,1]
3573              trA=(A11+A22)/2.
3574              A11-=trA
3575              A22-=trA
3576              s=sqrt(A12**2-A11*A22)
3577              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3578          elif s[0]==3:
3579              arg1=symmetric(arg)
3580              A11=arg1[0,0]
3581              A12=arg1[0,1]
3582              A22=arg1[1,1]
3583              A13=arg1[0,2]
3584              A23=arg1[1,2]
3585              A33=arg1[2,2]
3586              trA=(A11+A22+A33)/3.
3587              A11-=trA
3588              A22-=trA
3589              A33-=trA
3590              A13_2=A13**2
3591              A23_2=A23**2
3592              A12_2=A12**2
3593              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3594              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3595              sq_p=sqrt(p/3.)
3596              alpha_3=acos(clip(-q*(sq_p+whereZero(p,0.)*1.e-15)**(-3.)/2.,-1.,1.))/3.  # whereZero is protection against divison by zero
3597              sq_p*=2.
3598              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3599               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3600               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3601              return trA+sq_p*f
3602          else:
3603             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3604        elif isinstance(arg,float):
3605          return arg
3606        elif isinstance(arg,int):
3607          return float(arg)
3608        else:
3609          raise TypeError,"eigenvalues: Unknown argument type."
3610
3611    def eigenvalues_and_eigenvectors(arg):
3612        """
3613        returns the eigenvalues of the square matrix arg.
3614
3615        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3616                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3617        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3618        @return: the eigenvalues in increasing order.
3619        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3620        @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3621        """
3622        if isinstance(arg,numarray.NumArray):
3623          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3624        elif isinstance(arg,escript.Data):
3625          return arg._eigenvalues_and_eigenvectors()
3626        elif isinstance(arg,Symbol):
3627          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3628        elif isinstance(arg,float):
3629          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3630        elif isinstance(arg,int):
3631          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3632        else:
3633          raise TypeError,"eigenvalues: Unknown argument type."
3634  #=======================================================  #=======================================================
3635  #  Binary operations:  #  Binary operations:
3636  #=======================================================  #=======================================================
# Line 2798  def mult(arg0,arg1): Line 3749  def mult(arg0,arg1):
3749         """         """
3750         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3751         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3752            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3753         else:         else:
3754            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3755                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 2898  def quotient(arg0,arg1): Line 3849  def quotient(arg0,arg1):
3849         """         """
3850         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3851         if testForZero(args[0]):         if testForZero(args[0]):
3852            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3853         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3854            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3855               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3004  def power(arg0,arg1): Line 3955  def power(arg0,arg1):
3955         """         """
3956         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3957         if testForZero(args[0]):         if testForZero(args[0]):
3958            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3959         elif testForZero(args[1]):         elif testForZero(args[1]):
3960            return numarray.ones(args[0],numarray.Float)            return numarray.ones(pokeShape(args[1]),numarray.Float64)
3961         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3962            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3963         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):
# Line 3093  class Power_Symbol(DependendSymbol): Line 4044  class Power_Symbol(DependendSymbol):
4044           dargs=self.getDifferentiatedArguments(arg)           dargs=self.getDifferentiatedArguments(arg)
4046
4047  def maximum(arg0,arg1):  def maximum(*args):
4048      """      """
4049      the maximum of the two arguments      the maximum over arguments args
4050
4051      @param arg0: first argument      @param args: arguments
4052      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}      @type args: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4053      @param arg1: second argument      @return: is on object which gives at each entry the maximum of the coresponding values all args
@type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
@return: is on object which has the shape of arg0+arg1 and gives at each entry the maximum of the coresponding values of arg0 and arg1.
4054      @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
4055      """      """
4056      m=whereNegative(arg0-arg1)      out=None
4057      return m*arg1+(1.-m)*arg0      for a in args:
4058           if out==None:
4059              out=a
4060           else:
4063        return out
4064
4065  def minimum(arg0,arg1):  def minimum(*args):
4066      """      """
4067      the minumum of the two arguments      the minimum over arguments args
4068
4069      @param arg0: first argument      @param args: arguments
4070      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}      @type args: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4071      @param arg1: second argument      @return: is on object which gives at each entry the minimum of the coresponding values all args
@type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
@return: is on object which has the shape of arg0+arg1 and gives at each entry the minimum of the coresponding values of arg0 and arg1.
4072      @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
4073      """      """
4074      m=whereNegative(arg0-arg1)      out=None
4075      return m*arg0+(1.-m)*arg1      for a in args:
4076           if out==None:
4077              out=a
4078           else:
4081        return out
4082
4083    def clip(arg,minval=0.,maxval=1.):
4084        """
4085        cuts the values of arg between minval and maxval
4086
4087        @param arg: argument
4088        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4089        @param minval: lower range
4090        @type arg: C{float}
4091        @param maxval: upper range
4092        @type arg: C{float}
4093        @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4094                 less then maxval are unchanged.
4095        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4096        @raise ValueError: if minval>maxval
4097        """
4098        if minval>maxval:
4099           raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4100        return minimum(maximum(minval,arg),maxval)
4101
4102
4103  def inner(arg0,arg1):  def inner(arg0,arg1):
4104      """      """
# Line 3143  def inner(arg0,arg1): Line 4122  def inner(arg0,arg1):
4122      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4123      if not sh0==sh1:      if not sh0==sh1:
4124          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4125      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4126
4127  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4128      """      """
# Line 3171  def matrixmult(arg0,arg1): Line 4150  def matrixmult(arg0,arg1):
4150          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4151      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4152          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4153      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4154
4155  def outer(arg0,arg1):  def outer(arg0,arg1):
4156      """      """
# Line 3189  def outer(arg0,arg1): Line 4168  def outer(arg0,arg1):
4168      @return: the outer product of arg0 and arg1 at each data point      @return: the outer product of arg0 and arg1 at each data point
4169      @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
4170      """      """
4171      return generalTensorProduct(arg0,arg1,offset=0)      return generalTensorProduct(arg0,arg1,axis_offset=0)
4172
4173
4174  def tensormult(arg0,arg1):  def tensormult(arg0,arg1):
# Line 3231  def tensormult(arg0,arg1): Line 4210  def tensormult(arg0,arg1):
4210      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
4211      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4212      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4213         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4214      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):
4215         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4216      else:      else:
4217          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensormult: first argument must have rank 2 or 4"
4218
4219  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4220      """      """
4221      generalized tensor product      generalized tensor product
4222
4223      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4224
4225      where s runs through arg0.Shape[:arg0.Rank-offset]      where s runs through arg0.Shape[:arg0.Rank-axis_offset]
4226            r runs trough arg0.Shape[:offset]            r runs trough arg0.Shape[:axis_offset]
4227            t runs through arg1.Shape[offset:]            t runs through arg1.Shape[axis_offset:]
4228
4229      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
4230      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 3262  def generalTensorProduct(arg0,arg1,offse Line 4241  def generalTensorProduct(arg0,arg1,offse
4241      # 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
4242      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4243         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4244             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4245         else:         else:
4246             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4247                 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)
4248               arg0_c=arg0.copy()
4249               arg1_c=arg1.copy()
4250             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4251             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4252             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4253             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4254             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4255             arg0.resize((d0,d01))             arg0_c.resize((d0,d01))
4256             arg1.resize((d01,d1))             arg1_c.resize((d01,d1))
4257             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4258             for i0 in range(d0):             for i0 in range(d0):
4259                   for i1 in range(d1):                      for i1 in range(d1):
4260                        out[i0,i1]=numarray.sum(arg0[i0,:]*arg1[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4261             arg0.resize(sh0)             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
arg1.resize(sh1)
out.resize(sh0[:arg0.rank-offset]+sh1[offset:])
4262             return out             return out
4263      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4264         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4265             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4266         else:         else:
4267             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)
4268      else:            else:
4269         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4270
4271  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4272     """     """
4273     Symbol representing the quotient of two arguments.     Symbol representing the quotient of two arguments.
4274     """     """
4275     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4276         """         """
4277         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the quotient of two arguments
4278
# Line 3306  class GeneralTensorProduct_Symbol(Depend Line 4285  class GeneralTensorProduct_Symbol(Depend
4285         """         """
4286         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
4287         sh_arg1=pokeShape(arg1)         sh_arg1=pokeShape(arg1)
4288         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4289         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4290         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4291         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4292         if not sh01==sh10:         if not sh01==sh10:
4293             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)
4294         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])
4295
4296     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4297        """        """
# Line 3327  class GeneralTensorProduct_Symbol(Depend Line 4306  class GeneralTensorProduct_Symbol(Depend
4306        @raise: NotImplementedError: if the requested format is not available        @raise: NotImplementedError: if the requested format is not available
4307        """        """
4308        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4309           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])
4310        else:        else:
4311           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)
4312
# Line 3352  class GeneralTensorProduct_Symbol(Depend Line 4331  class GeneralTensorProduct_Symbol(Depend
4331           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4332           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4333
4334  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct
4335      "arg0 and arg1 are both Data objects but not neccesrily on the same function space!!!"      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4336      # calculate the return shape:      # calculate the return shape:
4337      shape0=arg0.getShape()[:arg0.getRank()-offset]      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4338      shape01=arg0.getShape()[arg0.getRank()-offset:]      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4339      shape10=arg1.getShape()[:offset]      shape10=arg1.getShape()[:axis_offset]
4340      shape1=arg1.getShape()[offset:]      shape1=arg1.getShape()[axis_offset:]
4341      if not shape01==shape10:      if not shape01==shape10:
4342          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)
4343
4344        # whatr function space should be used? (this here is not good!)
4345        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4346      # create return value:      # create return value:
4347      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())      out=escript.Data(0.,tuple(shape0+shape1),fs)
4348      #      #
4349      s0=[[]]      s0=[[]]
4350      for k in shape0:      for k in shape0:
# Line 3386  def escript_generalTensorProduct(arg0,ar Line 4367  def escript_generalTensorProduct(arg0,ar
4367
4368      for i0 in s0:      for i0 in s0:
4369         for i1 in s1:         for i1 in s1:
4370           s=escript.Scalar(0.,arg0.getFunctionSpace())           s=escript.Scalar(0.,fs)
4371           for i01 in s01:           for i01 in s01:
4372              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))
4373           out.__setitem__(tuple(i0+i1),s)           out.__setitem__(tuple(i0+i1),s)
4374      return out      return out
4375
4376
4377  #=========================================================  #=========================================================
4378  #   some little helpers  #  functions dealing with spatial dependency
4379  #=========================================================  #=========================================================
4381      """      """
4382      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4383
4384        If C{g} is the returned object, then
4385
4386      @param arg:   Data object representing the function which gradient        - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4387                    to be calculated.        - 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.
4388          - 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.
4389          - 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.
4390
4391        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4392        @type arg: L{escript.Data} or L{Symbol}
4393      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4394                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4395        @type where: C{None} or L{escript.FunctionSpace}
4397        @rtype:  L{escript.Data} or L{Symbol}
4398      """      """
4399      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4404         else:         else:
4406      else:      else:
4408
4410       """
4411       L{Symbol} representing the result of the gradient operator
4412       """
4413       def __init__(self,arg,where=None):
4414          """
4415          initialization of gradient L{Symbol} with argument arg
4416          @param arg: argument of function
4417          @type arg: L{Symbol}.
4418          @param where: FunctionSpace in which the gradient will be calculated.
4419                      If not present or C{None} an appropriate default is used.
4420          @type where: C{None} or L{escript.FunctionSpace}
4421          """
4422          d=arg.getDim()
4423          if d==None:
4424             raise ValueError,"argument must have a spatial dimension"
4426
4427       def getMyCode(self,argstrs,format="escript"):
4428          """
4429          returns a program code that can be used to evaluate the symbol.
4430
4431          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4432          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4433          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4434          @type format: C{str}
4435          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4436          @rtype: C{str}
4437          @raise: NotImplementedError: if the requested format is not available
4438          """
4439          if format=="escript" or format=="str"  or format=="text":
4441          else:
4442             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4443
4444       def substitute(self,argvals):
4445          """
4446          assigns new values to symbols in the definition of the symbol.
4447          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4448
4449          @param argvals: new values assigned to symbols
4450          @type argvals: C{dict} with keywords of type L{Symbol}.
4451          @return: result of the substitution process. Operations are executed as much as possible.
4452          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4453          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4454          """
4455          if argvals.has_key(self):
4456             arg=argvals[self]
4457             if self.isAppropriateValue(arg):
4458                return arg
4459             else:
4460                raise TypeError,"%s: new value is not appropriate."%str(self)
4461          else:
4462             arg=self.getSubstitutedArguments(argvals)
4464
4465       def diff(self,arg):
4466          """
4467          differential of this object
4468
4469          @param arg: the derivative is calculated with respect to arg
4470          @type arg: L{escript.Symbol}
4471          @return: derivative with respect to C{arg}
4472          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4473          """
4474          if arg==self:
4475             return identity(self.getShape())
4476          else:
4478
4479  def integrate(arg,where=None):  def integrate(arg,where=None):
4480      """      """
4481      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}
4482      its domain.      before integration.
4483
4484      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4485        @type arg: L{escript.Data} or L{Symbol}
4486      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4487                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4488        @type where: C{None} or L{escript.FunctionSpace}
4489        @return: integral of arg.
4490        @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}
4491      """      """
4492      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):
4493         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4494      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4495         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
# Line 3449  def integrate(arg,where=None): Line 4500  def integrate(arg,where=None):
4500      else:      else:
4501        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4502
4503    class Integrate_Symbol(DependendSymbol):
4504       """
4505       L{Symbol} representing the result of the spatial integration operator
4506       """
4507       def __init__(self,arg,where=None):
4508          """
4509          initialization of integration L{Symbol} with argument arg
4510          @param arg: argument of the integration
4511          @type arg: L{Symbol}.
4512          @param where: FunctionSpace in which the integration will be calculated.
4513                      If not present or C{None} an appropriate default is used.
4514          @type where: C{None} or L{escript.FunctionSpace}
4515          """
4516          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4517
4518       def getMyCode(self,argstrs,format="escript"):
4519          """
4520          returns a program code that can be used to evaluate the symbol.
4521
4522          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4523          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4524          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4525          @type format: C{str}
4526          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4527          @rtype: C{str}
4528          @raise: NotImplementedError: if the requested format is not available
4529          """
4530          if format=="escript" or format=="str"  or format=="text":
4531             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4532          else:
4533             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4534
4535       def substitute(self,argvals):
4536          """
4537          assigns new values to symbols in the definition of the symbol.
4538          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4539
4540          @param argvals: new values assigned to symbols
4541          @type argvals: C{dict} with keywords of type L{Symbol}.
4542          @return: result of the substitution process. Operations are executed as much as possible.
4543          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4544          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4545          """
4546          if argvals.has_key(self):
4547             arg=argvals[self]
4548             if self.isAppropriateValue(arg):
4549                return arg
4550             else:
4551                raise TypeError,"%s: new value is not appropriate."%str(self)
4552          else:
4553             arg=self.getSubstitutedArguments(argvals)
4554             return integrate(arg[0],where=arg[1])
4555
4556       def diff(self,arg):
4557          """
4558          differential of this object
4559
4560          @param arg: the derivative is calculated with respect to arg
4561          @type arg: L{escript.Symbol}
4562          @return: derivative with respect to C{arg}
4563          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4564          """
4565          if arg==self:
4566             return identity(self.getShape())
4567          else:
4568             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4569
4570
4571  def interpolate(arg,where):  def interpolate(arg,where):
4572      """      """
4573      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
4574
4575      @param arg:    interpolant      @param arg: interpolant
4576      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
4577        @param where: FunctionSpace to be interpolated to
4578        @type where: L{escript.FunctionSpace}
4579        @return: interpolated argument
4580        @rtype:  C{escript.Data} or L{Symbol}
4581      """      """
4582      if testForZero(arg):      if isinstance(arg,Symbol):
4583        return 0         return Interpolate_Symbol(arg,where)
elif isinstance(arg,Symbol):
return Interpolated_Symbol(arg,where)
4584      else:      else:
4585         return escript.Data(arg,where)         return escript.Data(arg,where)
4586
4587  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
4588      """     """
4589      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
4590       """
4591      @param arg:   Data object representing the function which gradient to     def __init__(self,arg,where):
4592                    be calculated.        """
4593      @param where: FunctionSpace in which the gradient will be calculated.        initialization of interpolation L{Symbol} with argument arg
4594                    If not present or C{None} an appropriate default is used.        @param arg: argument of the interpolation
4595      """        @type arg: L{Symbol}.
4596      g=grad(arg,where)        @param where: FunctionSpace into which the argument is interpolated.
4597      return trace(g,axis0=g.getRank()-2,axis1=g.getRank()-1)        @type where: L{escript.FunctionSpace}
4598          """
4599  def jump(arg):        super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
"""
Returns the jump of arg across a continuity.
4600
4601      @param arg:   Data object representing the function which gradient     def getMyCode(self,argstrs,format="escript"):
4602                    to be calculated.        """
4603      """        returns a program code that can be used to evaluate the symbol.
d=arg.getDomain()
return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())
4604
4605  #=============================        @param argstrs: gives for each argument a string representing the argument for the evaluation.
4606  #        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4607  # wrapper for various functions: if the argument has attribute the function name        @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4608  # as an argument it calls the corresponding methods. Otherwise the corresponding        @type format: C{str}
4609  # numarray function is called.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4610          @rtype: C{str}
4611          @raise: NotImplementedError: if the requested format is not available
4612          """
4613          if format=="escript" or format=="str"  or format=="text":
4614             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
4615          else:
4616             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4617
4618  # functions involving the underlying Domain:     def substitute(self,argvals):
4619          """
4620          assigns new values to symbols in the definition of the symbol.
4621          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4622
4623          @param argvals: new values assigned to symbols
4624          @type argvals: C{dict} with keywords of type L{Symbol}.
4625          @return: result of the substitution process. Operations are executed as much as possible.
4626          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4627          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4628          """
4629          if argvals.has_key(self):
4630             arg=argvals[self]
4631             if self.isAppropriateValue(arg):
4632                return arg
4633             else:
4634                raise TypeError,"%s: new value is not appropriate."%str(self)
4635          else:
4636             arg=self.getSubstitutedArguments(argvals)
4637             return interpolate(arg[0],where=arg[1])
4638
4639  # functions returning Data objects:     def diff(self,arg):
4640          """
4641          differential of this object
4642
4643  def minval(arg):        @param arg: the derivative is calculated with respect to arg
4644      return arg._minval()        @type arg: L{escript.Symbol}
4645          @return: derivative with respect to C{arg}
4646          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
4647          """
4648          if arg==self:
4649             return identity(self.getShape())
4650          else:
4651             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4652
def maxval(arg):
return arg._maxval()
4653
4654  def transpose(arg,axis=None):  def div(arg,where=None):
4655      """      """
4656      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
4657
4658      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
4659        @type arg: L{escript.Data} or L{Symbol}
4660        @param where: FunctionSpace in which the divergence will be calculated.
4661                      If not present or C{None} an appropriate default is used.
4662        @type where: C{None} or L{escript.FunctionSpace}
4663        @return: divergence of arg.
4664        @rtype:  L{escript.Data} or L{Symbol}
4665      """      """
if axis==None:
r=0
if hasattr(arg,"getRank"): r=arg.getRank()
if hasattr(arg,"rank"): r=arg.rank
axis=r/2
4666      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4667         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
4668      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4669         # 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)
4670      else:      else:
4671         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
4672        if not arg.getShape()==(dim,):
4673          raise ValueError,"div: expected shape is (%s,)"%dim
4675
4676  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
4677      """      """
4678      Return      returns the jump of arg across the continuity of the domain
4679
4680      @param arg:      @param arg: argument
4681        @type arg: L{escript.Data} or L{Symbol}
4682        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
4683                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
4684        @type domain: C{None} or L{escript.Domain}
4685        @return: jump of arg
4686        @rtype:  L{escript.Data} or L{Symbol}
4687      """      """
4688      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
4689         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)
4690
4691  def length(arg):  def L2(arg):
4692      """      """
4693        returns the L2 norm of arg at where
4694      @param arg:
4695        @param arg: function which L2 to be calculated.
4696        @type arg: L{escript.Data} or L{Symbol}
4697        @return: L2 norm of arg.
4698        @rtype:  L{float} or L{Symbol}
4699        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
4700      """      """
4701      if isinstance(arg,escript.Data):      return sqrt(integrate(inner(arg,arg)))
4702         if arg.isEmpty(): return escript.Data()  #=============================
4703         if arg.getRank()==0:  #
return abs(arg)
elif arg.getRank()==1:
out=escript.Scalar(0,arg.getFunctionSpace())
for i in range(arg.getShape()[0]):
out+=arg[i]**2
return sqrt(out)
elif arg.getRank()==2:
out=escript.Scalar(0,arg.getFunctionSpace())
for i in range(arg.getShape()[0]):
for j in range(arg.getShape()[1]):
out+=arg[i,j]**2
return sqrt(out)
elif arg.getRank()==3:
out=escript.Scalar(0,arg.getFunctionSpace())
for i in range(arg.getShape()[0]):
for j in range(arg.getShape()[1]):
for k in range(arg.getShape()[2]):
out+=arg[i,j,k]**2
return sqrt(out)
elif arg.getRank()==4:
out=escript.Scalar(0,arg.getFunctionSpace())
for i in range(arg.getShape()[0]):
for j in range(arg.getShape()[1]):
for k in range(arg.getShape()[2]):
for l in range(arg.getShape()[3]):
out+=arg[i,j,k,l]**2
return sqrt(out)
else:
raise SystemError,"length is not been fully implemented yet"
# return arg.length()
elif isinstance(arg,float):
return abs(arg)
else:
return sqrt((arg**2).sum())
4704
4705  def reorderComponents(arg,index):  def reorderComponents(arg,index):
4706      """      """
4707      resorts the component of arg according to index      resorts the component of arg according to index
4708
4709      """      """
4710      pass      raise NotImplementedError
4711  #  #
4712  # \$Log: util.py,v \$  # \$Log: util.py,v \$
4713  # 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.290 changed lines Added in v.587