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

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

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

revision 612 by gross, Tue Mar 21 23:54:39 2006 UTC revision 1042 by gross, Mon Mar 19 03:50:34 2007 UTC
# Line 3  Line 3 
3  """  """
4  Utility functions for escript  Utility functions for escript
5    
 @remark:  This module is under construction and is still tested!!!  
   
6  @var __author__: name of author  @var __author__: name of author
7  @var __licence__: licence agreement  @var __copyright__: copyrights
8    @var __license__: licence agreement
9  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
10  @var __version__: version  @var __version__: version
11  @var __date__: date of the version  @var __date__: date of the version
# Line 16  __author__="Lutz Gross, l.gross@uq.edu.a Line 15  __author__="Lutz Gross, l.gross@uq.edu.a
15  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                      http://www.access.edu.au                      http://www.access.edu.au
17                  Primary Business: Queensland, Australia"""                  Primary Business: Queensland, Australia"""
18  __licence__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19               http://www.opensource.org/licences/osl-3.0.php"""               http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
21  __version__="$Revision$"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
# Line 27  import math Line 26  import math
26  import numarray  import numarray
27  import escript  import escript
28  import os  import os
29    from esys.escript import C_GeneralTensorProduct
30    
31  #=========================================================  #=========================================================
32  #   some helpers:  #   some helpers:
33  #=========================================================  #=========================================================
34    def getTagNames(domain):
35        """
36        returns a list of the tag names used by the domain
37    
38        
39        @param domain: a domain object
40        @type domain: C{escript.Domain}
41        @return: a list of the tag name used by the domain.
42        @rtype: C{list} of C{str}
43        """
44        return [n.strip() for n in domain.showTagNames().split(",") ]
45        
46  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
47      """      """
48      writes a L{Data} objects into a files using the the VTK XML file format.      writes a L{Data} objects into a files using the the VTK XML file format.
49    
50      Example:      Example::
51    
52         tmp=Scalar(..)         tmp=Scalar(..)
53         v=Vector(..)         v=Vector(..)
# Line 63  def saveDX(filename,domain=None,**data): Line 75  def saveDX(filename,domain=None,**data):
75      """      """
76      writes a L{Data} objects into a files using the the DX file format.      writes a L{Data} objects into a files using the the DX file format.
77    
78      Example:      Example::
79    
80         tmp=Scalar(..)         tmp=Scalar(..)
81         v=Vector(..)         v=Vector(..)
# Line 94  def kronecker(d=3): Line 106  def kronecker(d=3):
106     @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
107     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
108     @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
109     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2.
110     """     """
111     return identityTensor(d)     return identityTensor(d)
112    
# Line 130  def identityTensor(d=3): Line 142  def identityTensor(d=3):
142     @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
143     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
144     @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
145     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
146     """     """
147     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
148         return escript.Data(identity((d.getDim(),)),d)         return escript.Data(identity((d.getDim(),)),d)
# Line 146  def identityTensor4(d=3): Line 158  def identityTensor4(d=3):
158     @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
159     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
160     @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
161     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
162     """     """
163     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
164         return escript.Data(identity((d.getDim(),d.getDim())),d)         return escript.Data(identity((d.getDim(),d.getDim())),d)
# Line 164  def unitVector(i=0,d=3): Line 176  def unitVector(i=0,d=3):
176     @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
177     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
178     @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
179     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 1     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
180     """     """
181     return kronecker(d)[i]     return kronecker(d)[i]
182    
# Line 220  def inf(arg): Line 232  def inf(arg):
232    
233      @param arg: argument      @param arg: argument
234      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
235      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
236      @rtype: C{float}      @rtype: C{float}
237      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
238      """      """
# Line 239  def inf(arg): Line 251  def inf(arg):
251  #=========================================================================  #=========================================================================
252  #   some little helpers  #   some little helpers
253  #=========================================================================  #=========================================================================
254  def pokeShape(arg):  def getRank(arg):
255        """
256        identifies the rank of its argument
257    
258        @param arg: a given object
259        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
260        @return: the rank of the argument
261        @rtype: C{int}
262        @raise TypeError: if type of arg cannot be processed
263        """
264    
265        if isinstance(arg,numarray.NumArray):
266            return arg.rank
267        elif isinstance(arg,escript.Data):
268            return arg.getRank()
269        elif isinstance(arg,float):
270            return 0
271        elif isinstance(arg,int):
272            return 0
273        elif isinstance(arg,Symbol):
274            return arg.getRank()
275        else:
276          raise TypeError,"getShape: cannot identify shape"
277    def getShape(arg):
278      """      """
279      identifies the shape of its argument      identifies the shape of its argument
280    
# Line 261  def pokeShape(arg): Line 296  def pokeShape(arg):
296      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
297          return arg.getShape()          return arg.getShape()
298      else:      else:
299        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
300    
301  def pokeDim(arg):  def pokeDim(arg):
302      """      """
# Line 284  def commonShape(arg0,arg1): Line 319  def commonShape(arg0,arg1):
319      """      """
320      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.
321    
322      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
323      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
324      @return: the shape of arg0 or arg1 such that the left port equals the shape of arg0 and the right end equals the shape of arg1.      @return: the shape of arg0 or arg1 such that the left port equals the shape of arg0 and the right end equals the shape of arg1.
325      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
326      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
327      """      """
328      sh0=pokeShape(arg0)      sh0=getShape(arg0)
329      sh1=pokeShape(arg1)      sh1=getShape(arg1)
330      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
331         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
332               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
# Line 309  def commonDim(*args): Line 344  def commonDim(*args):
344      """      """
345      identifies, if possible, the spatial dimension across a set of objects which may or my not have a spatial dimension.      identifies, if possible, the spatial dimension across a set of objects which may or my not have a spatial dimension.
346    
347      @param *args: given objects      @param args: given objects
348      @return: the spatial dimension of the objects with identifiable dimension (see L{pokeDim}). If none the objects has      @return: the spatial dimension of the objects with identifiable dimension (see L{pokeDim}). If none the objects has
349               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
350      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 331  def testForZero(arg): Line 366  def testForZero(arg):
366    
367      @param arg: a given object      @param arg: a given object
368      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}
369      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
370      @rtype : C{bool}      @rtype: C{bool}
371      """      """
372      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
373         return not Lsup(arg)>0.         return not Lsup(arg)>0.
# Line 435  def matchType(arg0=0.,arg1=0.): Line 470  def matchType(arg0=0.,arg1=0.):
470    
471  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
472      """      """
473            return representations of arg0 amd arg1 which ahve the same shape
474    
475      If shape is not given the shape "largest" shape of args is used.      @param arg0: a given object
476        @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
477      @param args: a given ob      @param arg1: a given object
478      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
479      @return: True if the argument is identical to zero.      @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
480      @rtype: C{list} of C{int}      @rtype: C{tuple}
481      """      """
482      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
483      sh0=pokeShape(arg0)      sh0=getShape(arg0)
484      sh1=pokeShape(arg1)      sh1=getShape(arg1)
485      if len(sh0)<len(sh):      if len(sh0)<len(sh):
486         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
487      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
# Line 470  class Symbol(object): Line 505  class Symbol(object):
505         Creates an instance of a symbol of a given shape. The symbol may depending on a list of arguments args which may be         Creates an instance of a symbol of a given shape. The symbol may depending on a list of arguments args which may be
506         symbols or any other object.         symbols or any other object.
507    
508         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
509         @type arg: C{list}         @type args: C{list}
510         @param shape: the shape         @param shape: the shape
511         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
512         @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
# Line 514  class Symbol(object): Line 549  class Symbol(object):
549         """         """
550         the shape of the symbol.         the shape of the symbol.
551    
552         @return : the shape of the symbol.         @return: the shape of the symbol.
553         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
554         """         """
555         return self.__shape         return self.__shape
# Line 523  class Symbol(object): Line 558  class Symbol(object):
558         """         """
559         the spatial dimension         the spatial dimension
560    
561         @return : the spatial dimension         @return: the spatial dimension
562         @rtype: C{int} if the dimension is defined. Otherwise C{None} is returned.         @rtype: C{int} if the dimension is defined. Otherwise C{None} is returned.
563         """         """
564         return self.__dim         return self.__dim
# Line 547  class Symbol(object): Line 582  class Symbol(object):
582         """         """
583         substitutes symbols in the arguments of this object and returns the result as a list.         substitutes symbols in the arguments of this object and returns the result as a list.
584    
585         @param argvals: L{Symbols} and their substitutes. The L{Symbol} u in the expression defining this object is replaced by argvals[u].         @param argvals: L{Symbol} and their substitutes. The L{Symbol} u in the expression defining this object is replaced by argvals[u].
586         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
587         @rtype: C{list} of objects         @rtype: C{list} of objects
588         @return: list of the object assigned to the arguments through substitution or for the arguments which are not L{Symbols} the value assigned to the argument at instantiation.         @return: list of the object assigned to the arguments through substitution or for the arguments which are not L{Symbol} the value assigned to the argument at instantiation.
589         """         """
590         out=[]         out=[]
591         for a in self.getArgument():         for a in self.getArgument():
# Line 574  class Symbol(object): Line 609  class Symbol(object):
609            if isinstance(a,Symbol):            if isinstance(a,Symbol):
610               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
611            else:            else:
612                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
613                if len(s)>0:                if len(s)>0:
614                   out.append(numarray.zeros(s),numarray.Float64)                   out.append(numarray.zeros(s),numarray.Float64)
615                else:                else:
# Line 674  class Symbol(object): Line 709  class Symbol(object):
709         """         """
710         returns -self.         returns -self.
711    
712         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
713         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
714         """         """
715         return self*(-1.)         return self*(-1.)
# Line 683  class Symbol(object): Line 718  class Symbol(object):
718         """         """
719         returns +self.         returns +self.
720    
721         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
722         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
723         """         """
724         return self*(1.)         return self*(1.)
725    
726     def __abs__(self):     def __abs__(self):
727         """         """
728         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
729         """         """
730         return Abs_Symbol(self)         return Abs_Symbol(self)
731    
# Line 700  class Symbol(object): Line 735  class Symbol(object):
735    
736         @param other: object to be added to this object         @param other: object to be added to this object
737         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
738         @return:  a S{Symbol} representing the sum of this object and C{other}         @return:  a L{Symbol} representing the sum of this object and C{other}
739         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
740         """         """
741         return add(self,other)         return add(self,other)
# Line 711  class Symbol(object): Line 746  class Symbol(object):
746    
747         @param other: object this object is added to         @param other: object this object is added to
748         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
749         @return: a S{Symbol} representing the sum of C{other} and this object object         @return: a L{Symbol} representing the sum of C{other} and this object object
750         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
751         """         """
752         return add(other,self)         return add(other,self)
# Line 722  class Symbol(object): Line 757  class Symbol(object):
757    
758         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
759         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
760         @return: a S{Symbol} representing the difference of C{other} and this object         @return: a L{Symbol} representing the difference of C{other} and this object
761         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
762         """         """
763         return add(self,-other)         return add(self,-other)
# Line 733  class Symbol(object): Line 768  class Symbol(object):
768    
769         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
770         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
771         @return: a S{Symbol} representing the difference of this object and C{other}.         @return: a L{Symbol} representing the difference of this object and C{other}.
772         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
773         """         """
774         return add(-self,other)         return add(-self,other)
# Line 744  class Symbol(object): Line 779  class Symbol(object):
779    
780         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
781         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
782         @return: a S{Symbol} representing the product of the object and C{other}.         @return: a L{Symbol} representing the product of the object and C{other}.
783         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
784         """         """
785         return mult(self,other)         return mult(self,other)
# Line 755  class Symbol(object): Line 790  class Symbol(object):
790    
791         @param other: object this object is multiplied with         @param other: object this object is multiplied with
792         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
793         @return: a S{Symbol} representing the product of C{other} and the object.         @return: a L{Symbol} representing the product of C{other} and the object.
794         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
795         """         """
796         return mult(other,self)         return mult(other,self)
# Line 766  class Symbol(object): Line 801  class Symbol(object):
801    
802         @param other: object dividing this object         @param other: object dividing this object
803         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
804         @return: a S{Symbol} representing the quotient of this object and C{other}         @return: a L{Symbol} representing the quotient of this object and C{other}
805         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
806         """         """
807         return quotient(self,other)         return quotient(self,other)
# Line 777  class Symbol(object): Line 812  class Symbol(object):
812    
813         @param other: object dividing this object         @param other: object dividing this object
814         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
815         @return: a S{Symbol} representing the quotient of C{other} and this object         @return: a L{Symbol} representing the quotient of C{other} and this object
816         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
817         """         """
818         return quotient(other,self)         return quotient(other,self)
# Line 788  class Symbol(object): Line 823  class Symbol(object):
823    
824         @param other: exponent         @param other: exponent
825         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
826         @return: a S{Symbol} representing the power of this object to C{other}         @return: a L{Symbol} representing the power of this object to C{other}
827         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
828         """         """
829         return power(self,other)         return power(self,other)
# Line 799  class Symbol(object): Line 834  class Symbol(object):
834    
835         @param other: basis         @param other: basis
836         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
837         @return: a S{Symbol} representing the power of C{other} to this object         @return: a L{Symbol} representing the power of C{other} to this object
838         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
839         """         """
840         return power(other,self)         return power(other,self)
# Line 810  class Symbol(object): Line 845  class Symbol(object):
845    
846         @param index: defines a         @param index: defines a
847         @type index: C{slice} or C{int} or a C{tuple} of them         @type index: C{slice} or C{int} or a C{tuple} of them
848         @return: a S{Symbol} representing the slice defined by index         @return: a L{Symbol} representing the slice defined by index
849         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
850         """         """
851         return GetSlice_Symbol(self,index)         return GetSlice_Symbol(self,index)
# Line 820  class DependendSymbol(Symbol): Line 855  class DependendSymbol(Symbol):
855     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.
856     Two DependendSymbol are equal if they have the same shape, the same arguments and one of them has an unspecified spatial dimension or the spatial dimension is identical       Two DependendSymbol are equal if they have the same shape, the same arguments and one of them has an unspecified spatial dimension or the spatial dimension is identical  
857        
858     Example:     Example::
859        
860     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
861     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
862     print u1==u2       print u1==u2
863     False       False
864        
865        but       but::
866    
867     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
868     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
869     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
870     print u1==u2, u1==u3       print u1==u2, u1==u3
871     True False       True False
872    
873     @note: DependendSymbol should be used as return value of functions with L{Symbol} arguments. This will allow the optimizer to remove redundant function calls.     @note: DependendSymbol should be used as return value of functions with L{Symbol} arguments. This will allow the optimizer to remove redundant function calls.
874     """     """
# Line 923  class GetSlice_Symbol(DependendSymbol): Line 958  class GetSlice_Symbol(DependendSymbol):
958        @type format: C{str}        @type format: C{str}
959        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
960        @rtype: C{str}        @rtype: C{str}
961        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
962        """        """
963        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
964           return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])           return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
# Line 959  def log10(arg): Line 994  def log10(arg):
994    
995     @param arg: argument     @param arg: argument
996     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
997     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
998     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
999     """     """
1000     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 981  def wherePositive(arg): Line 1016  def wherePositive(arg):
1016    
1017     @param arg: argument     @param arg: argument
1018     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1019     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1020     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1021     """     """
1022     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1027  class WherePositive_Symbol(DependendSymb Line 1062  class WherePositive_Symbol(DependendSymb
1062        @type format: C{str}        @type format: C{str}
1063        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1064        @rtype: C{str}        @rtype: C{str}
1065        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1066        """        """
1067        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1068            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1063  def whereNegative(arg): Line 1098  def whereNegative(arg):
1098    
1099     @param arg: argument     @param arg: argument
1100     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1101     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1102     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1103     """     """
1104     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1109  class WhereNegative_Symbol(DependendSymb Line 1144  class WhereNegative_Symbol(DependendSymb
1144        @type format: C{str}        @type format: C{str}
1145        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1146        @rtype: C{str}        @rtype: C{str}
1147        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1148        """        """
1149        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1150            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1145  def whereNonNegative(arg): Line 1180  def whereNonNegative(arg):
1180    
1181     @param arg: argument     @param arg: argument
1182     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1183     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1184     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1185     """     """
1186     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1175  def whereNonPositive(arg): Line 1210  def whereNonPositive(arg):
1210    
1211     @param arg: argument     @param arg: argument
1212     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1213     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1214     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1215     """     """
1216     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1207  def whereZero(arg,tol=0.): Line 1242  def whereZero(arg,tol=0.):
1242     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1243     @param tol: tolerance. values with absolute value less then tol are accepted as zero.     @param tol: tolerance. values with absolute value less then tol are accepted as zero.
1244     @type tol: C{float}     @type tol: C{float}
1245     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1246     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1247     """     """
1248     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1215  def whereZero(arg,tol=0.): Line 1250  def whereZero(arg,tol=0.):
1250        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1251        return out        return out
1252     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1253        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1254     elif isinstance(arg,float):     elif isinstance(arg,float):
1255        if abs(arg)<=tol:        if abs(arg)<=tol:
1256          return 1.          return 1.
# Line 1256  class WhereZero_Symbol(DependendSymbol): Line 1288  class WhereZero_Symbol(DependendSymbol):
1288        @type format: C{str}        @type format: C{str}
1289        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1290        @rtype: C{str}        @rtype: C{str}
1291        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1292        """        """
1293        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1294           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1290  def whereNonZero(arg,tol=0.): Line 1322  def whereNonZero(arg,tol=0.):
1322    
1323     @param arg: argument     @param arg: argument
1324     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1325     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1326     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1327     """     """
1328     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1298  def whereNonZero(arg,tol=0.): Line 1330  def whereNonZero(arg,tol=0.):
1330        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1331        return out        return out
1332     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1333        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1334     elif isinstance(arg,float):     elif isinstance(arg,float):
1335        if abs(arg)>tol:        if abs(arg)>tol:
1336          return 1.          return 1.
# Line 1317  def whereNonZero(arg,tol=0.): Line 1346  def whereNonZero(arg,tol=0.):
1346     else:     else:
1347        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1348    
1349    def erf(arg):
1350       """
1351       returns erf of argument arg
1352    
1353       @param arg: argument
1354       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1355       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1356       @raises TypeError: if the type of the argument is not expected.
1357       """
1358       if isinstance(arg,escript.Data):
1359          return arg._erf()
1360       else:
1361          raise TypeError,"erf: Unknown argument type."
1362    
1363  def sin(arg):  def sin(arg):
1364     """     """
1365     returns sine of argument arg     returns sine of argument arg
1366    
1367     @param arg: argument     @param arg: argument
1368     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1369     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1370     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1371     """     """
1372     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1361  class Sin_Symbol(DependendSymbol): Line 1404  class Sin_Symbol(DependendSymbol):
1404        @type format: C{str}        @type format: C{str}
1405        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1406        @rtype: C{str}        @rtype: C{str}
1407        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1408        """        """
1409        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1410            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1413  def cos(arg): Line 1456  def cos(arg):
1456    
1457     @param arg: argument     @param arg: argument
1458     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1459     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1460     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1461     """     """
1462     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1451  class Cos_Symbol(DependendSymbol): Line 1494  class Cos_Symbol(DependendSymbol):
1494        @type format: C{str}        @type format: C{str}
1495        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1496        @rtype: C{str}        @rtype: C{str}
1497        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1498        """        """
1499        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1500            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1503  def tan(arg): Line 1546  def tan(arg):
1546    
1547     @param arg: argument     @param arg: argument
1548     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1549     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1550     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1551     """     """
1552     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1541  class Tan_Symbol(DependendSymbol): Line 1584  class Tan_Symbol(DependendSymbol):
1584        @type format: C{str}        @type format: C{str}
1585        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1586        @rtype: C{str}        @rtype: C{str}
1587        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1588        """        """
1589        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1590            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1593  def asin(arg): Line 1636  def asin(arg):
1636    
1637     @param arg: argument     @param arg: argument
1638     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1639     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1640     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1641     """     """
1642     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1631  class Asin_Symbol(DependendSymbol): Line 1674  class Asin_Symbol(DependendSymbol):
1674        @type format: C{str}        @type format: C{str}
1675        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1676        @rtype: C{str}        @rtype: C{str}
1677        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1678        """        """
1679        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1680            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1683  def acos(arg): Line 1726  def acos(arg):
1726    
1727     @param arg: argument     @param arg: argument
1728     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1729     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1730     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1731     """     """
1732     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1721  class Acos_Symbol(DependendSymbol): Line 1764  class Acos_Symbol(DependendSymbol):
1764        @type format: C{str}        @type format: C{str}
1765        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1766        @rtype: C{str}        @rtype: C{str}
1767        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1768        """        """
1769        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1770            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1773  def atan(arg): Line 1816  def atan(arg):
1816    
1817     @param arg: argument     @param arg: argument
1818     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1819     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1820     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1821     """     """
1822     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1811  class Atan_Symbol(DependendSymbol): Line 1854  class Atan_Symbol(DependendSymbol):
1854        @type format: C{str}        @type format: C{str}
1855        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1856        @rtype: C{str}        @rtype: C{str}
1857        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1858        """        """
1859        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1860            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1863  def sinh(arg): Line 1906  def sinh(arg):
1906    
1907     @param arg: argument     @param arg: argument
1908     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1909     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1910     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1911     """     """
1912     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1901  class Sinh_Symbol(DependendSymbol): Line 1944  class Sinh_Symbol(DependendSymbol):
1944        @type format: C{str}        @type format: C{str}
1945        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1946        @rtype: C{str}        @rtype: C{str}
1947        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1948        """        """
1949        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1950            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1953  def cosh(arg): Line 1996  def cosh(arg):
1996    
1997     @param arg: argument     @param arg: argument
1998     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1999     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2000     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2001     """     """
2002     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1991  class Cosh_Symbol(DependendSymbol): Line 2034  class Cosh_Symbol(DependendSymbol):
2034        @type format: C{str}        @type format: C{str}
2035        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2036        @rtype: C{str}        @rtype: C{str}
2037        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2038        """        """
2039        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2040            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2043  def tanh(arg): Line 2086  def tanh(arg):
2086    
2087     @param arg: argument     @param arg: argument
2088     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2089     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2090     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2091     """     """
2092     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2081  class Tanh_Symbol(DependendSymbol): Line 2124  class Tanh_Symbol(DependendSymbol):
2124        @type format: C{str}        @type format: C{str}
2125        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2126        @rtype: C{str}        @rtype: C{str}
2127        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2128        """        """
2129        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2130            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2133  def asinh(arg): Line 2176  def asinh(arg):
2176    
2177     @param arg: argument     @param arg: argument
2178     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2179     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2180     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2181     """     """
2182     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2171  class Asinh_Symbol(DependendSymbol): Line 2214  class Asinh_Symbol(DependendSymbol):
2214        @type format: C{str}        @type format: C{str}
2215        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2216        @rtype: C{str}        @rtype: C{str}
2217        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2218        """        """
2219        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2220            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2223  def acosh(arg): Line 2266  def acosh(arg):
2266    
2267     @param arg: argument     @param arg: argument
2268     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2269     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2270     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2271     """     """
2272     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2261  class Acosh_Symbol(DependendSymbol): Line 2304  class Acosh_Symbol(DependendSymbol):
2304        @type format: C{str}        @type format: C{str}
2305        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2306        @rtype: C{str}        @rtype: C{str}
2307        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2308        """        """
2309        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2310            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2313  def atanh(arg): Line 2356  def atanh(arg):
2356    
2357     @param arg: argument     @param arg: argument
2358     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2359     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2360     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2361     """     """
2362     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2351  class Atanh_Symbol(DependendSymbol): Line 2394  class Atanh_Symbol(DependendSymbol):
2394        @type format: C{str}        @type format: C{str}
2395        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2396        @rtype: C{str}        @rtype: C{str}
2397        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2398        """        """
2399        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2400            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2403  def exp(arg): Line 2446  def exp(arg):
2446    
2447     @param arg: argument     @param arg: argument
2448     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2449     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2450     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2451     """     """
2452     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2441  class Exp_Symbol(DependendSymbol): Line 2484  class Exp_Symbol(DependendSymbol):
2484        @type format: C{str}        @type format: C{str}
2485        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2486        @rtype: C{str}        @rtype: C{str}
2487        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2488        """        """
2489        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2490            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2493  def sqrt(arg): Line 2536  def sqrt(arg):
2536    
2537     @param arg: argument     @param arg: argument
2538     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2539     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2540     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2541     """     """
2542     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2531  class Sqrt_Symbol(DependendSymbol): Line 2574  class Sqrt_Symbol(DependendSymbol):
2574        @type format: C{str}        @type format: C{str}
2575        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2576        @rtype: C{str}        @rtype: C{str}
2577        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2578        """        """
2579        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2580            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2583  def log(arg): Line 2626  def log(arg):
2626    
2627     @param arg: argument     @param arg: argument
2628     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2629     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2630     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2631     """     """
2632     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2621  class Log_Symbol(DependendSymbol): Line 2664  class Log_Symbol(DependendSymbol):
2664        @type format: C{str}        @type format: C{str}
2665        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2666        @rtype: C{str}        @rtype: C{str}
2667        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2668        """        """
2669        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2670            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2673  def sign(arg): Line 2716  def sign(arg):
2716    
2717     @param arg: argument     @param arg: argument
2718     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2719     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2720     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2721     """     """
2722     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2721  class Abs_Symbol(DependendSymbol): Line 2764  class Abs_Symbol(DependendSymbol):
2764        @type format: C{str}        @type format: C{str}
2765        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2766        @rtype: C{str}        @rtype: C{str}
2767        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2768        """        """
2769        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2770            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2773  def minval(arg): Line 2816  def minval(arg):
2816    
2817     @param arg: argument     @param arg: argument
2818     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2819     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2820     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2821     """     """
2822     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2814  class Minval_Symbol(DependendSymbol): Line 2857  class Minval_Symbol(DependendSymbol):
2857        @type format: C{str}        @type format: C{str}
2858        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2859        @rtype: C{str}        @rtype: C{str}
2860        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2861        """        """
2862        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2863            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2850  def maxval(arg): Line 2893  def maxval(arg):
2893    
2894     @param arg: argument     @param arg: argument
2895     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2896     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2897     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2898     """     """
2899     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2891  class Maxval_Symbol(DependendSymbol): Line 2934  class Maxval_Symbol(DependendSymbol):
2934        @type format: C{str}        @type format: C{str}
2935        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2936        @rtype: C{str}        @rtype: C{str}
2937        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2938        """        """
2939        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2940            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2927  def length(arg): Line 2970  def length(arg):
2970    
2971     @param arg: argument     @param arg: argument
2972     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2973     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2974     """     """
2975     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2976    
# Line 2937  def trace(arg,axis_offset=0): Line 2980  def trace(arg,axis_offset=0):
2980    
2981     @param arg: argument     @param arg: argument
2982     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2983     @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     @param axis_offset: C{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
2984                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
2985     @type axis_offset: C{int}     @type axis_offset: C{int}
2986     @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.     @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2987     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
# Line 2946  def trace(arg,axis_offset=0): Line 2989  def trace(arg,axis_offset=0):
2989     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
2990        sh=arg.shape        sh=arg.shape
2991        if len(sh)<2:        if len(sh)<2:
2992          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2993        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
2994          raise ValueError,"trace: axis_offset must be between 0 and %s"%len(sh)-2          raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
2995        s1=1        s1=1
2996        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
2997        s2=1        s2=1
2998        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2999        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
3000          raise ValueError,"trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3001        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
3002        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
3003        for i1 in range(s1):        for i1 in range(s1):
# Line 2963  def trace(arg,axis_offset=0): Line 3006  def trace(arg,axis_offset=0):
3006        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3007        return out        return out
3008     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3009        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
3010            raise ValueError,"rank of argument must be greater than 1"
3011          if axis_offset<0 or axis_offset>arg.getRank()-2:
3012            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3013          s=list(arg.getShape())        
3014          if not s[axis_offset] == s[axis_offset+1]:
3015            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3016          return arg._trace(axis_offset)
3017     elif isinstance(arg,float):     elif isinstance(arg,float):
3018        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
3019     elif isinstance(arg,int):     elif isinstance(arg,int):
3020        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
3021     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3022        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
3023     else:     else:
3024        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
3025    
 def escript_trace(arg,axis_offset): # this should be escript._trace  
       "arg si a Data objects!!!"  
       if arg.getRank()<2:  
         raise ValueError,"escript_trace: rank of argument must be greater than 1"  
       if axis_offset<0 or axis_offset>arg.getRank()-2:  
         raise ValueError,"escript_trace: axis_offset must be between 0 and %s"%arg.getRank()-2  
       s=list(arg.getShape())          
       if not s[axis_offset] == s[axis_offset+1]:  
         raise ValueError,"escript_trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)  
       out=escript.Data(0.,tuple(s[0:axis_offset]+s[axis_offset+2:]),arg.getFunctionSpace())  
       if arg.getRank()==2:  
          for i0 in range(s[0]):  
             out+=arg[i0,i0]  
       elif arg.getRank()==3:  
          if axis_offset==0:  
             for i0 in range(s[0]):  
                   for i2 in range(s[2]):  
                          out[i2]+=arg[i0,i0,i2]  
          elif axis_offset==1:  
             for i0 in range(s[0]):  
                for i1 in range(s[1]):  
                          out[i0]+=arg[i0,i1,i1]  
       elif arg.getRank()==4:  
          if axis_offset==0:  
             for i0 in range(s[0]):  
                   for i2 in range(s[2]):  
                      for i3 in range(s[3]):  
                          out[i2,i3]+=arg[i0,i0,i2,i3]  
          elif axis_offset==1:  
             for i0 in range(s[0]):  
                for i1 in range(s[1]):  
                      for i3 in range(s[3]):  
                          out[i0,i3]+=arg[i0,i1,i1,i3]  
          elif axis_offset==2:  
             for i0 in range(s[0]):  
                for i1 in range(s[1]):  
                   for i2 in range(s[2]):  
                          out[i0,i1]+=arg[i0,i1,i2,i2]  
       return out  
3026  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
3027     """     """
3028     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 3021  class Trace_Symbol(DependendSymbol): Line 3032  class Trace_Symbol(DependendSymbol):
3032        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3033        @param arg: argument of function        @param arg: argument of function
3034        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3035        @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        @param axis_offset: C{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
3036                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3037        @type axis_offset: C{int}        @type axis_offset: C{int}
3038        """        """
3039        if arg.getRank()<2:        if arg.getRank()<2:
3040          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3041        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3042          raise ValueError,"Trace_Symbol: axis_offset must be between 0 and %s"%arg.getRank()-2          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3043        s=list(arg.getShape())                s=list(arg.getShape())        
3044        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3045          raise ValueError,"Trace_Symbol: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3046        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3047    
3048     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 3044  class Trace_Symbol(DependendSymbol): Line 3055  class Trace_Symbol(DependendSymbol):
3055        @type format: C{str}        @type format: C{str}
3056        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3057        @rtype: C{str}        @rtype: C{str}
3058        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3059        """        """
3060        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3061           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3088  class Trace_Symbol(DependendSymbol): Line 3099  class Trace_Symbol(DependendSymbol):
3099    
3100  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3101     """     """
3102     returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.     returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3103    
3104     @param arg: argument     @param arg: argument
3105     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3106     @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.     @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3107                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3108     @type axis_offset: C{int}     @type axis_offset: C{int}
3109     @return: transpose of arg     @return: transpose of arg
3110     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
# Line 3102  def transpose(arg,axis_offset=None): Line 3113  def transpose(arg,axis_offset=None):
3113        if axis_offset==None: axis_offset=int(arg.rank/2)        if axis_offset==None: axis_offset=int(arg.rank/2)
3114        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3115     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3116        if axis_offset==None: axis_offset=int(arg.getRank()/2)        r=arg.getRank()
3117        return escript_transpose(arg,axis_offset)        if axis_offset==None: axis_offset=int(r/2)
3118          if axis_offset<0 or axis_offset>r:
3119            raise ValueError,"axis_offset must be between 0 and %s"%r
3120          return arg._transpose(axis_offset)
3121     elif isinstance(arg,float):     elif isinstance(arg,float):
3122        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3123          raise ValueError,"transpose: axis_offset must be 0 for float argument"          raise ValueError,"axis_offset must be 0 for float argument"
3124        return arg        return arg
3125     elif isinstance(arg,int):     elif isinstance(arg,int):
3126        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3127          raise ValueError,"transpose: axis_offset must be 0 for int argument"          raise ValueError,"axis_offset must be 0 for int argument"
3128        return float(arg)        return float(arg)
3129     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3130        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3131        return Transpose_Symbol(arg,axis_offset)        return Transpose_Symbol(arg,axis_offset)
3132     else:     else:
3133        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"Unknown argument type."
3134    
 def escript_transpose(arg,axis_offset): # this should be escript._transpose  
       "arg si a Data objects!!!"  
       r=arg.getRank()  
       if axis_offset<0 or axis_offset>r:  
         raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r  
       s=arg.getShape()  
       s_out=s[axis_offset:]+s[:axis_offset]  
       out=escript.Data(0.,s_out,arg.getFunctionSpace())  
       if r==4:  
          if axis_offset==1:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i3,i0,i1,i2]  
          elif axis_offset==2:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i2,i3,i0,i1]  
          elif axis_offset==3:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i1,i2,i3,i0]  
          else:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i0,i1,i2,i3]  
       elif r==3:  
          if axis_offset==1:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                          out[i0,i1,i2]=arg[i2,i0,i1]  
          elif axis_offset==2:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                          out[i0,i1,i2]=arg[i1,i2,i0]  
          else:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                          out[i0,i1,i2]=arg[i0,i1,i2]  
       elif r==2:  
          if axis_offset==1:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                          out[i0,i1]=arg[i1,i0]  
          else:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                          out[i0,i1]=arg[i0,i1]  
       elif r==1:  
           for i0 in range(s_out[0]):  
                out[i0]=arg[i0]  
       elif r==0:  
              out=arg+0.  
       return out  
3135  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3136     """     """
3137     L{Symbol} representing the result of the transpose function     L{Symbol} representing the result of the transpose function
# Line 3192  class Transpose_Symbol(DependendSymbol): Line 3142  class Transpose_Symbol(DependendSymbol):
3142    
3143        @param arg: argument of function        @param arg: argument of function
3144        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3145         @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.        @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3146                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3147        @type axis_offset: C{int}        @type axis_offset: C{int}
3148        """        """
3149        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3150        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3151          raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3152        s=arg.getShape()        s=arg.getShape()
3153        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3154    
# Line 3212  class Transpose_Symbol(DependendSymbol): Line 3162  class Transpose_Symbol(DependendSymbol):
3162        @type format: C{str}        @type format: C{str}
3163        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3164        @rtype: C{str}        @rtype: C{str}
3165        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3166        """        """
3167        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3168           return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3253  class Transpose_Symbol(DependendSymbol): Line 3203  class Transpose_Symbol(DependendSymbol):
3203           return identity(self.getShape())           return identity(self.getShape())
3204        else:        else:
3205           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3206    
3207    def swap_axes(arg,axis0=0,axis1=1):
3208       """
3209       returns the swap of arg by swaping the components axis0 and axis1
3210    
3211       @param arg: argument
3212       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3213       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3214       @type axis0: C{int}
3215       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3216       @type axis1: C{int}
3217       @return: C{arg} with swaped components
3218       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3219       """
3220       if axis0 > axis1:
3221          axis0,axis1=axis1,axis0
3222       if isinstance(arg,numarray.NumArray):
3223          return numarray.swapaxes(arg,axis0,axis1)
3224       elif isinstance(arg,escript.Data):
3225          return arg._swap_axes(axis0,axis1)
3226       elif isinstance(arg,float):
3227          raise TyepError,"float argument is not supported."
3228       elif isinstance(arg,int):
3229          raise TyepError,"int argument is not supported."
3230       elif isinstance(arg,Symbol):
3231          return SwapAxes_Symbol(arg,axis0,axis1)
3232       else:
3233          raise TypeError,"Unknown argument type."
3234    
3235    class SwapAxes_Symbol(DependendSymbol):
3236       """
3237       L{Symbol} representing the result of the swap function
3238       """
3239       def __init__(self,arg,axis0=0,axis1=1):
3240          """
3241          initialization of swap L{Symbol} with argument arg
3242    
3243          @param arg: argument
3244          @type arg: L{Symbol}.
3245          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3246          @type axis0: C{int}
3247          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3248          @type axis1: C{int}
3249          """
3250          if arg.getRank()<2:
3251             raise ValueError,"argument must have at least rank 2."
3252          if axis0<0 or axis0>arg.getRank()-1:
3253             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3254          if axis1<0 or axis1>arg.getRank()-1:
3255             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3256          if axis0 == axis1:
3257             raise ValueError,"axis indices must be different."
3258          if axis0 > axis1:
3259             axis0,axis1=axis1,axis0
3260          s=arg.getShape()
3261          s_out=[]
3262          for i in range(len(s)):
3263             if i == axis0:
3264                s_out.append(s[axis1])
3265             elif i == axis1:
3266                s_out.append(s[axis0])
3267             else:
3268                s_out.append(s[i])
3269          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3270    
3271       def getMyCode(self,argstrs,format="escript"):
3272          """
3273          returns a program code that can be used to evaluate the symbol.
3274    
3275          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3276          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3277          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3278          @type format: C{str}
3279          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3280          @rtype: C{str}
3281          @raise NotImplementedError: if the requested format is not available
3282          """
3283          if format=="escript" or format=="str"  or format=="text":
3284             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3285          else:
3286             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3287    
3288       def substitute(self,argvals):
3289          """
3290          assigns new values to symbols in the definition of the symbol.
3291          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3292    
3293          @param argvals: new values assigned to symbols
3294          @type argvals: C{dict} with keywords of type L{Symbol}.
3295          @return: result of the substitution process. Operations are executed as much as possible.
3296          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3297          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3298          """
3299          if argvals.has_key(self):
3300             arg=argvals[self]
3301             if self.isAppropriateValue(arg):
3302                return arg
3303             else:
3304                raise TypeError,"%s: new value is not appropriate."%str(self)
3305          else:
3306             arg=self.getSubstitutedArguments(argvals)
3307             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3308    
3309       def diff(self,arg):
3310          """
3311          differential of this object
3312    
3313          @param arg: the derivative is calculated with respect to arg
3314          @type arg: L{escript.Symbol}
3315          @return: derivative with respect to C{arg}
3316          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3317          """
3318          if arg==self:
3319             return identity(self.getShape())
3320          else:
3321             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3322    
3323  def symmetric(arg):  def symmetric(arg):
3324      """      """
3325      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
# Line 3265  def symmetric(arg): Line 3332  def symmetric(arg):
3332      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3333        if arg.rank==2:        if arg.rank==2:
3334          if not (arg.shape[0]==arg.shape[1]):          if not (arg.shape[0]==arg.shape[1]):
3335             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3336        elif arg.rank==4:        elif arg.rank==4:
3337          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3338             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3339        else:        else:
3340          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3341        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3342      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3343        return escript_symmetric(arg)        if arg.getRank()==2:
3344            if not (arg.getShape()[0]==arg.getShape()[1]):
3345               raise ValueError,"argument must be square."
3346            return arg._symmetric()
3347          elif arg.getRank()==4:
3348            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3349               raise ValueError,"argument must be square."
3350            return arg._symmetric()
3351          else:
3352            raise ValueError,"rank 2 or 4 is required."
3353      elif isinstance(arg,float):      elif isinstance(arg,float):
3354        return arg        return arg
3355      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3281  def symmetric(arg): Line 3357  def symmetric(arg):
3357      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3358        if arg.getRank()==2:        if arg.getRank()==2:
3359          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3360             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3361        elif arg.getRank()==4:        elif arg.getRank()==4:
3362          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3363             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3364        else:        else:
3365          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3366        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3367      else:      else:
3368        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3369    
 def escript_symmetric(arg): # this should be implemented in c++  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               out[i0,i1]=(arg[i0,i1]+arg[i1,i0])/2.  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               for i2 in range(arg.getShape()[2]):  
                  for i3 in range(arg.getShape()[3]):  
                      out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]+arg[i2,i3,i0,i1])/2.  
       else:  
         raise ValueError,"escript_symmetric: rank 2 or 4 is required."  
       return out  
   
3370  def nonsymmetric(arg):  def nonsymmetric(arg):
3371      """      """
3372      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
# Line 3332  def nonsymmetric(arg): Line 3387  def nonsymmetric(arg):
3387          raise ValueError,"nonsymmetric: rank 2 or 4 is required."          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3388        return (arg-transpose(arg))/2        return (arg-transpose(arg))/2
3389      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3390        return escript_nonsymmetric(arg)        if arg.getRank()==2:
3391            if not (arg.getShape()[0]==arg.getShape()[1]):
3392               raise ValueError,"argument must be square."
3393            return arg._nonsymmetric()
3394          elif arg.getRank()==4:
3395            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3396               raise ValueError,"argument must be square."
3397            return arg._nonsymmetric()
3398          else:
3399            raise ValueError,"rank 2 or 4 is required."
3400      elif isinstance(arg,float):      elif isinstance(arg,float):
3401        return arg        return arg
3402      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3350  def nonsymmetric(arg): Line 3414  def nonsymmetric(arg):
3414      else:      else:
3415        raise TypeError,"nonsymmetric: Unknown argument type."        raise TypeError,"nonsymmetric: Unknown argument type."
3416    
 def escript_nonsymmetric(arg): # this should be implemented in c++  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               out[i0,i1]=(arg[i0,i1]-arg[i1,i0])/2.  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               for i2 in range(arg.getShape()[2]):  
                  for i3 in range(arg.getShape()[3]):  
                      out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]-arg[i2,i3,i0,i1])/2.  
       else:  
         raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."  
       return out  
   
   
3417  def inverse(arg):  def inverse(arg):
3418      """      """
3419      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3420    
3421      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3422      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3423      @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])      @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3424      @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
3425      @remark: for L{escript.Data} objects the dimension is restricted to 3.      @note: for L{escript.Data} objects the dimension is restricted to 3.
3426      """      """
3427      import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.      import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3428      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
# Line 3475  class Inverse_Symbol(DependendSymbol): Line 3517  class Inverse_Symbol(DependendSymbol):
3517        @type format: C{str}        @type format: C{str}
3518        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3519        @rtype: C{str}        @rtype: C{str}
3520        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3521        """        """
3522        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3523           return "inverse(%s)"%argstrs[0]           return "inverse(%s)"%argstrs[0]
# Line 3515  class Inverse_Symbol(DependendSymbol): Line 3557  class Inverse_Symbol(DependendSymbol):
3557        if arg==self:        if arg==self:
3558           return identity(self.getShape())           return identity(self.getShape())
3559        else:        else:
3560           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3561    
3562  def eigenvalues(arg):  def eigenvalues(arg):
3563      """      """
# Line 3526  def eigenvalues(arg): Line 3568  def eigenvalues(arg):
3568      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3569      @return: the eigenvalues in increasing order.      @return: the eigenvalues in increasing order.
3570      @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.
3571      @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.      @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3572      """      """
3573      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3574        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
# Line 3596  def eigenvalues_and_eigenvectors(arg): Line 3638  def eigenvalues_and_eigenvectors(arg):
3638               eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is               eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3639               the eigenvector coresponding to the i-th eigenvalue.               the eigenvector coresponding to the i-th eigenvalue.
3640      @rtype: L{tuple} of L{escript.Data}.      @rtype: L{tuple} of L{escript.Data}.
3641      @remark: The dimension is restricted to 3.      @note: The dimension is restricted to 3.
3642      """      """
3643      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3644        raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"        raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
# Line 3653  class Add_Symbol(DependendSymbol): Line 3695  class Add_Symbol(DependendSymbol):
3695         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3696         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3697         """         """
3698         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3699         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3700         if not sh0==sh1:         if not sh0==sh1:
3701            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3702         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3669  class Add_Symbol(DependendSymbol): Line 3711  class Add_Symbol(DependendSymbol):
3711        @type format: C{str}        @type format: C{str}
3712        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3713        @rtype: C{str}        @rtype: C{str}
3714        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3715        """        """
3716        if format=="str" or format=="text":        if format=="str" or format=="text":
3717           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 3728  def mult(arg0,arg1): Line 3770  def mult(arg0,arg1):
3770         """         """
3771         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3772         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3773            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3774         else:         else:
3775            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3776                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3752  class Mult_Symbol(DependendSymbol): Line 3794  class Mult_Symbol(DependendSymbol):
3794         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3795         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3796         """         """
3797         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3798         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3799         if not sh0==sh1:         if not sh0==sh1:
3800            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3801         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3768  class Mult_Symbol(DependendSymbol): Line 3810  class Mult_Symbol(DependendSymbol):
3810        @type format: C{str}        @type format: C{str}
3811        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3812        @rtype: C{str}        @rtype: C{str}
3813        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3814        """        """
3815        if format=="str" or format=="text":        if format=="str" or format=="text":
3816           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3828  def quotient(arg0,arg1): Line 3870  def quotient(arg0,arg1):
3870         """         """
3871         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3872         if testForZero(args[0]):         if testForZero(args[0]):
3873            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3874         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3875            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3876               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3857  class Quotient_Symbol(DependendSymbol): Line 3899  class Quotient_Symbol(DependendSymbol):
3899         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3900         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3901         """         """
3902         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3903         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3904         if not sh0==sh1:         if not sh0==sh1:
3905            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3906         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3873  class Quotient_Symbol(DependendSymbol): Line 3915  class Quotient_Symbol(DependendSymbol):
3915        @type format: C{str}        @type format: C{str}
3916        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3917        @rtype: C{str}        @rtype: C{str}
3918        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3919        """        """
3920        if format=="str" or format=="text":        if format=="str" or format=="text":
3921           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3934  def power(arg0,arg1): Line 3976  def power(arg0,arg1):
3976         """         """
3977         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3978         if testForZero(args[0]):         if testForZero(args[0]):
3979            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3980         elif testForZero(args[1]):         elif testForZero(args[1]):
3981            return numarray.ones(pokeShape(args[1]),numarray.Float64)            return numarray.ones(getShape(args[1]),numarray.Float64)
3982         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3983            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3984         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 3959  class Power_Symbol(DependendSymbol): Line 4001  class Power_Symbol(DependendSymbol):
4001         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4002         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4003         """         """
4004         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4005         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4006         if not sh0==sh1:         if not sh0==sh1:
4007            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4008         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3977  class Power_Symbol(DependendSymbol): Line 4019  class Power_Symbol(DependendSymbol):
4019        @type format: C{str}        @type format: C{str}
4020        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4021        @rtype: C{str}        @rtype: C{str}
4022        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4023        """        """
4024        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4025           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 4059  def minimum(*args): Line 4101  def minimum(*args):
4101            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4102      return out      return out
4103    
4104  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4105      """      """
4106      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4107    
4108      @param arg: argument      @param arg: argument
4109      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4110      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4111      @type arg: C{float}      @type minval: C{float} or C{None}
4112      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4113      @type arg: C{float}      @type maxval: C{float} or C{None}
4114      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4115               less then maxval are unchanged.               less then maxval are unchanged.
4116      @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
4117      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4118      """      """
4119      if minval>maxval:      if not minval==None and not maxval==None:
4120         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4121      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4122        if minval == None:
4123            tmp=arg
4124        else:
4125            tmp=maximum(minval,arg)
4126        if maxval == None:
4127            return tmp
4128        else:
4129            return minimum(tmp,maxval)
4130    
4131        
4132  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 4093  def inner(arg0,arg1): Line 4143  def inner(arg0,arg1):
4143      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4144      @param arg1: second argument      @param arg1: second argument
4145      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4146      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4147      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input
4148      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4149      """      """
4150      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4151      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4152      if not sh0==sh1:      if not sh0==sh1:
4153          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4154      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4155    
4156    def outer(arg0,arg1):
4157        """
4158        the outer product of the two argument:
4159    
4160        out[t,s]=arg0[t]*arg1[s]
4161    
4162        where
4163    
4164            - s runs through arg0.Shape
4165            - t runs through arg1.Shape
4166    
4167        @param arg0: first argument
4168        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4169        @param arg1: second argument
4170        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4171        @return: the outer product of arg0 and arg1 at each data point
4172        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4173        """
4174        return generalTensorProduct(arg0,arg1,axis_offset=0)
4175    
4176  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4177      """      """
4178        see L{matrix_mult}
4179        """
4180        return matrix_mult(arg0,arg1)
4181    
4182    def matrix_mult(arg0,arg1):
4183        """
4184      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4185    
4186      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4187    
4188            or      or
4189    
4190      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4191    
4192      The second dimension of arg0 and the length of arg1 must match.      The second dimension of arg0 and the first dimension of arg1 must match.
4193    
4194      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4195      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4123  def matrixmult(arg0,arg1): Line 4199  def matrixmult(arg0,arg1):
4199      @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
4200      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4201      """      """
4202      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4203      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4204      if not len(sh0)==2 :      if not len(sh0)==2 :
4205          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4206      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4207          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4208      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4209    
4210  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4211      """      """
4212      the outer product of the two argument:      use L{tensor_mult}
   
     out[t,s]=arg0[t]*arg1[s]  
   
     where s runs through arg0.Shape  
           t runs through arg1.Shape  
   
     @param arg0: first argument  
     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}  
     @param arg1: second argument  
     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}  
     @return: the outer product of arg0 and arg1 at each data point  
     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input  
4213      """      """
4214      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
   
4215    
4216  def tensormult(arg0,arg1):  def tensor_mult(arg0,arg1):
4217      """      """
4218      the tensor product of the two argument:      the tensor product of the two argument:
   
4219            
4220      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4221    
4222      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4223    
4224                   or      or
4225    
4226      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4227    
# Line 4168  def tensormult(arg0,arg1): Line 4230  def tensormult(arg0,arg1):
4230    
4231      out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2,s3]      out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2,s3]
4232                                
4233                   or      or
4234    
4235      out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2]      out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2]
4236    
4237                   or      or
4238    
4239      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]
4240    
4241      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 last dimension of arg1 must match and  
4242      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 two first dimensions of arg1.
4243    
4244      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4245      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4186  def tensormult(arg0,arg1): Line 4248  def tensormult(arg0,arg1):
4248      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4249      @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
4250      """      """
4251      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4252      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4253      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4254         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4255      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):
4256         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4257      else:      else:
4258          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4259    
4260  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4261      """      """
# Line 4201  def generalTensorProduct(arg0,arg1,axis_ Line 4263  def generalTensorProduct(arg0,arg1,axis_
4263    
4264      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4265    
4266      where s runs through arg0.Shape[:arg0.Rank-axis_offset]      where
           r runs trough arg0.Shape[:axis_offset]  
           t runs through arg1.Shape[axis_offset:]  
4267    
4268      In the first case the the second dimension of arg0 and the length of arg1 must match and            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4269      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4270            - t runs through arg1.Shape[axis_offset:]
4271    
4272      @param arg0: first argument      @param arg0: first argument
4273      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4274      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4275      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4276      @return: the general tensor product of arg0 and arg1 at each data point.      @return: the general tensor product of arg0 and arg1 at each data point.
4277      @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
# Line 4223  def generalTensorProduct(arg0,arg1,axis_ Line 4284  def generalTensorProduct(arg0,arg1,axis_
4284             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4285         else:         else:
4286             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4287                 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)                 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)
4288             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4289             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4290             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 4249  def generalTensorProduct(arg0,arg1,axis_ Line 4310  def generalTensorProduct(arg0,arg1,axis_
4310                                    
4311  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4312     """     """
4313     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4314     """     """
4315     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4316         """         """
4317         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4318    
4319         @param arg0: numerator         @param arg0: first argument
4320         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4321         @param arg1: denominator         @param arg1: second argument
4322         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4323         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4324         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4325         """         """
4326         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4327         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4328         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4329         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4330         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 4282  class GeneralTensorProduct_Symbol(Depend Line 4343  class GeneralTensorProduct_Symbol(Depend
4343        @type format: C{str}        @type format: C{str}
4344        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4345        @rtype: C{str}        @rtype: C{str}
4346        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4347        """        """
4348        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4349           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
# Line 4310  class GeneralTensorProduct_Symbol(Depend Line 4371  class GeneralTensorProduct_Symbol(Depend
4371           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4372           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4373    
4374  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4375      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4376      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4377      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
4378      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  def transposed_matrix_mult(arg0,arg1):
4379      shape10=arg1.getShape()[:axis_offset]      """
4380      shape1=arg1.getShape()[axis_offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4381      if not shape01==shape10:  
4382          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)      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
4383    
4384      # whatr function space should be used? (this here is not good!)      or
4385      fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
4386      # create return value:      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4387      out=escript.Data(0.,tuple(shape0+shape1),fs)  
4388      #      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4389      s0=[[]]  
4390      for k in shape0:      The first dimension of arg0 and arg1 must match.
4391            s=[]  
4392            for j in s0:      @param arg0: first argument of rank 2
4393                  for i in range(k): s.append(j+[slice(i,i)])      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4394            s0=s      @param arg1: second argument of at least rank 1
4395      s1=[[]]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4396      for k in shape1:      @return: the product of the transposed of arg0 and arg1 at each data point
4397            s=[]      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4398            for j in s1:      @raise ValueError: if the shapes of the arguments are not appropriate
4399                  for i in range(k): s.append(j+[slice(i,i)])      """
4400            s1=s      sh0=getShape(arg0)
4401      s01=[[]]      sh1=getShape(arg1)
4402      for k in shape01:      if not len(sh0)==2 :
4403            s=[]          raise ValueError,"first argument must have rank 2"
4404            for j in s01:      if not len(sh1)==2 and not len(sh1)==1:
4405                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"second argument must have rank 1 or 2"
4406            s01=s      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4407    
4408      for i0 in s0:  def transposed_tensor_mult(arg0,arg1):
4409         for i1 in s1:      """
4410           s=escript.Scalar(0.,fs)      the tensor product of the transposed of the first and the second argument
4411           for i01 in s01:      
4412              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      for arg0 of rank 2 this is
4413           out.__setitem__(tuple(i0+i1),s)  
4414      return out      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4415    
4416        or
4417    
4418        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4419    
4420      
4421        and for arg0 of rank 4 this is
4422    
4423        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4424                  
4425        or
4426    
4427        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4428    
4429        or
4430    
4431        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4432    
4433        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4434        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4435    
4436        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4437    
4438        @param arg0: first argument of rank 2 or 4
4439        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4440        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4441        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4442        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4443        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4444        """
4445        sh0=getShape(arg0)
4446        sh1=getShape(arg1)
4447        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4448           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4449        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4450           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4451        else:
4452            raise ValueError,"first argument must have rank 2 or 4"
4453    
4454    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4455        """
4456        generalized tensor product of transposed of arg0 and arg1:
4457    
4458        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4459    
4460        where
4461    
4462            - s runs through arg0.Shape[axis_offset:]
4463            - r runs trough arg0.Shape[:axis_offset]
4464            - t runs through arg1.Shape[axis_offset:]
4465    
4466        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4467        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4468    
4469        @param arg0: first argument
4470        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4471        @param arg1: second argument
4472        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4473        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4474        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4475        """
4476        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4477        arg0,arg1=matchType(arg0,arg1)
4478        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4479        if isinstance(arg0,numarray.NumArray):
4480           if isinstance(arg1,Symbol):
4481               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4482           else:
4483               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4484                   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)
4485               arg0_c=arg0.copy()
4486               arg1_c=arg1.copy()
4487               sh0,sh1=arg0.shape,arg1.shape
4488               d0,d1,d01=1,1,1
4489               for i in sh0[axis_offset:]: d0*=i
4490               for i in sh1[axis_offset:]: d1*=i
4491               for i in sh0[:axis_offset]: d01*=i
4492               arg0_c.resize((d01,d0))
4493               arg1_c.resize((d01,d1))
4494               out=numarray.zeros((d0,d1),numarray.Float64)
4495               for i0 in range(d0):
4496                        for i1 in range(d1):
4497                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4498               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4499               return out
4500        elif isinstance(arg0,escript.Data):
4501           if isinstance(arg1,Symbol):
4502               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4503           else:
4504               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4505        else:      
4506           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4507                    
4508    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4509       """
4510       Symbol representing the general tensor product of the transposed of arg0 and arg1
4511       """
4512       def __init__(self,arg0,arg1,axis_offset=0):
4513           """
4514           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4515    
4516           @param arg0: first argument
4517           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4518           @param arg1: second argument
4519           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4520           @raise ValueError: inconsistent dimensions of arguments.
4521           @note: if both arguments have a spatial dimension, they must equal.
4522           """
4523           sh_arg0=getShape(arg0)
4524           sh_arg1=getShape(arg1)
4525           sh01=sh_arg0[:axis_offset]
4526           sh10=sh_arg1[:axis_offset]
4527           sh0=sh_arg0[axis_offset:]
4528           sh1=sh_arg1[axis_offset:]
4529           if not sh01==sh10:
4530               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)
4531           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4532    
4533       def getMyCode(self,argstrs,format="escript"):
4534          """
4535          returns a program code that can be used to evaluate the symbol.
4536    
4537          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4538          @type argstrs: C{list} of length 2 of C{str}.
4539          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4540          @type format: C{str}
4541          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4542          @rtype: C{str}
4543          @raise NotImplementedError: if the requested format is not available
4544          """
4545          if format=="escript" or format=="str" or format=="text":
4546             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4547          else:
4548             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4549    
4550       def substitute(self,argvals):
4551          """
4552          assigns new values to symbols in the definition of the symbol.
4553          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4554    
4555          @param argvals: new values assigned to symbols
4556          @type argvals: C{dict} with keywords of type L{Symbol}.
4557          @return: result of the substitution process. Operations are executed as much as possible.
4558          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4559          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4560          """
4561          if argvals.has_key(self):
4562             arg=argvals[self]
4563             if self.isAppropriateValue(arg):
4564                return arg
4565             else:
4566                raise TypeError,"%s: new value is not appropriate."%str(self)
4567          else:
4568             args=self.getSubstitutedArguments(argvals)
4569             return generalTransposedTensorProduct(args[0],args[1],args[2])
4570    
4571    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4572        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4573        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4574    
4575    def matrix_transposed_mult(arg0,arg1):
4576        """
4577        matrix-transposed(matrix) product of the two argument:
4578    
4579        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4580    
4581        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4582    
4583        The last dimensions of arg0 and arg1 must match.
4584    
4585        @param arg0: first argument of rank 2
4586        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4587        @param arg1: second argument of rank 2
4588        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4589        @return: the product of arg0 and the transposed of arg1 at each data point
4590        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4591        @raise ValueError: if the shapes of the arguments are not appropriate
4592        """
4593        sh0=getShape(arg0)
4594        sh1=getShape(arg1)
4595        if not len(sh0)==2 :
4596            raise ValueError,"first argument must have rank 2"
4597        if not len(sh1)==2 and not len(sh1)==1:
4598            raise ValueError,"second argument must have rank 1 or 2"
4599        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4600    
4601    def tensor_transposed_mult(arg0,arg1):
4602        """
4603        the tensor product of the first and the transpose of the second argument
4604        
4605        for arg0 of rank 2 this is
4606    
4607        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4608    
4609        and for arg0 of rank 4 this is
4610    
4611        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4612                  
4613        or
4614    
4615        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4616    
4617        In the first case the the second dimension of arg0 and arg1 must match and  
4618        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4619    
4620        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4621    
4622        @param arg0: first argument of rank 2 or 4
4623        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4624        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4625        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4626        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4627        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4628        """
4629        sh0=getShape(arg0)
4630        sh1=getShape(arg1)
4631        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4632           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4633        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4634           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4635        else:
4636            raise ValueError,"first argument must have rank 2 or 4"
4637    
4638    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4639        """
4640        generalized tensor product of transposed of arg0 and arg1:
4641    
4642        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4643    
4644        where
4645    
4646            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4647            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4648            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4649    
4650        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4651        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4652    
4653        @param arg0: first argument
4654        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4655        @param arg1: second argument
4656        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4657        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4658        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4659        """
4660        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4661        arg0,arg1=matchType(arg0,arg1)
4662        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4663        if isinstance(arg0,numarray.NumArray):
4664           if isinstance(arg1,Symbol):
4665               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4666           else:
4667               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4668                   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)
4669               arg0_c=arg0.copy()
4670               arg1_c=arg1.copy()
4671               sh0,sh1=arg0.shape,arg1.shape
4672               d0,d1,d01=1,1,1
4673               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4674               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4675               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4676               arg0_c.resize((d0,d01))
4677               arg1_c.resize((d1,d01))
4678               out=numarray.zeros((d0,d1),numarray.Float64)
4679               for i0 in range(d0):
4680                        for i1 in range(d1):
4681                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4682               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4683               return out
4684        elif isinstance(arg0,escript.Data):
4685           if isinstance(arg1,Symbol):
4686               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4687           else:
4688               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4689        else:      
4690           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4691                    
4692    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4693       """
4694       Symbol representing the general tensor product of arg0 and the transpose of arg1
4695       """
4696       def __init__(self,arg0,arg1,axis_offset=0):
4697           """
4698           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4699    
4700           @param arg0: first argument
4701           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4702           @param arg1: second argument
4703           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4704           @raise ValueError: inconsistent dimensions of arguments.
4705           @note: if both arguments have a spatial dimension, they must equal.
4706           """
4707           sh_arg0=getShape(arg0)
4708           sh_arg1=getShape(arg1)
4709           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4710           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4711           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4712           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4713           if not sh01==sh10:
4714               raise ValueError,"dimensions of last %s components in left argument don't match the last %s components in the right argument."%(axis_offset,axis_offset)
4715           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4716    
4717       def getMyCode(self,argstrs,format="escript"):
4718          """
4719          returns a program code that can be used to evaluate the symbol.
4720    
4721          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4722          @type argstrs: C{list} of length 2 of C{str}.
4723          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4724          @type format: C{str}
4725          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4726          @rtype: C{str}
4727          @raise NotImplementedError: if the requested format is not available
4728          """
4729          if format=="escript" or format=="str" or format=="text":
4730             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4731          else:
4732             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4733    
4734       def substitute(self,argvals):
4735          """
4736          assigns new values to symbols in the definition of the symbol.
4737          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4738    
4739          @param argvals: new values assigned to symbols
4740          @type argvals: C{dict} with keywords of type L{Symbol}.
4741          @return: result of the substitution process. Operations are executed as much as possible.
4742          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4743          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4744          """
4745          if argvals.has_key(self):
4746             arg=argvals[self]
4747             if self.isAppropriateValue(arg):
4748                return arg
4749             else:
4750                raise TypeError,"%s: new value is not appropriate."%str(self)
4751          else:
4752             args=self.getSubstitutedArguments(argvals)
4753             return generalTensorTransposedProduct(args[0],args[1],args[2])
4754    
4755    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4756        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4757        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4758    
4759  #=========================================================  #=========================================================
4760  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
# Line 4373  def grad(arg,where=None): Line 4776  def grad(arg,where=None):
4776                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4777      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4778      @return: gradient of arg.      @return: gradient of arg.
4779      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4780      """      """
4781      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4782         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 4413  class Grad_Symbol(DependendSymbol): Line 4816  class Grad_Symbol(DependendSymbol):
4816        @type format: C{str}        @type format: C{str}
4817        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4818        @rtype: C{str}        @rtype: C{str}
4819        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4820        """        """
4821        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4822           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4466  def integrate(arg,where=None): Line 4869  def integrate(arg,where=None):
4869                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4870      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4871      @return: integral of arg.      @return: integral of arg.
4872      @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}      @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4873      """      """
4874      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4875         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
# Line 4504  class Integrate_Symbol(DependendSymbol): Line 4907  class Integrate_Symbol(DependendSymbol):
4907        @type format: C{str}        @type format: C{str}
4908        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4909        @rtype: C{str}        @rtype: C{str}
4910        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4911        """        """
4912        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4913           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4556  def interpolate(arg,where): Line 4959  def interpolate(arg,where):
4959      @param where: FunctionSpace to be interpolated to      @param where: FunctionSpace to be interpolated to
4960      @type where: L{escript.FunctionSpace}      @type where: L{escript.FunctionSpace}
4961      @return: interpolated argument      @return: interpolated argument
4962      @rtype:  C{escript.Data} or L{Symbol}      @rtype: C{escript.Data} or L{Symbol}
4963      """      """
4964      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4965         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
# Line 4587  class Interpolate_Symbol(DependendSymbol Line 4990  class Interpolate_Symbol(DependendSymbol
4990        @type format: C{str}        @type format: C{str}
4991        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4992        @rtype: C{str}        @rtype: C{str}
4993        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4994        """        """
4995        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4996           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4640  def div(arg,where=None): Line 5043  def div(arg,where=None):
5043                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
5044      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
5045      @return: divergence of arg.      @return: divergence of arg.
5046      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5047      """      """
5048      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5049          dim=arg.getDim()          dim=arg.getDim()
# Line 4662  def jump(arg,domain=None): Line 5065  def jump(arg,domain=None):
5065                     the domain of arg is used. If arg is a L{Symbol} the domain must be present.                     the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5066      @type domain: C{None} or L{escript.Domain}      @type domain: C{None} or L{escript.Domain}
5067      @return: jump of arg      @return: jump of arg
5068      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5069      """      """
5070      if domain==None: domain=arg.getDomain()      if domain==None: domain=arg.getDomain()
5071      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
# Line 4674  def L2(arg): Line 5077  def L2(arg):
5077      @param arg: function which L2 to be calculated.      @param arg: function which L2 to be calculated.
5078      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
5079      @return: L2 norm of arg.      @return: L2 norm of arg.
5080      @rtype:  L{float} or L{Symbol}      @rtype: L{float} or L{Symbol}
5081      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5082      """      """
5083      return sqrt(integrate(inner(arg,arg)))      return sqrt(integrate(inner(arg,arg)))

Legend:
Removed from v.612  
changed lines
  Added in v.1042

  ViewVC Help
Powered by ViewVC 1.1.26