/[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 588 by gross, Fri Mar 10 04:45:04 2006 UTC revision 881 by gross, Thu Oct 26 02:54:47 2006 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2    
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
12  """  """
13                                                                                                                                                                                                                                                                                                                                                                                                            
14  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
15  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                        http://www.access.edu.au
17                    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19                 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 30  __date__="$Date$" Line 24  __date__="$Date$"
24    
25  import math  import math
26  import numarray  import numarray
 import numarray.linear_algebra  
27  import escript  import escript
28  import os  import os
29    from esys.escript import C_GeneralTensorProduct
 # missing tests:  
   
 # def pokeShape(arg):  
 # def pokeDim(arg):  
 # def commonShape(arg0,arg1):  
 # def commonDim(*args):  
 # def testForZero(arg):  
 # def matchType(arg0=0.,arg1=0.):  
 # def matchShape(arg0,arg1):  
   
 # def reorderComponents(arg,index):  
   
 #  
 # slicing: get  
 #          set  
 #  
 # and derivatives  
30    
31  #=========================================================  #=========================================================
32  #   some helpers:  #   some helpers:
# Line 59  def saveVTK(filename,domain=None,**data) Line 35  def saveVTK(filename,domain=None,**data)
35      """      """
36      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.
37    
38      Example:      Example::
39    
40         tmp=Scalar(..)         tmp=Scalar(..)
41         v=Vector(..)         v=Vector(..)
# Line 87  def saveDX(filename,domain=None,**data): Line 63  def saveDX(filename,domain=None,**data):
63      """      """
64      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.
65    
66      Example:      Example::
67    
68         tmp=Scalar(..)         tmp=Scalar(..)
69         v=Vector(..)         v=Vector(..)
# Line 118  def kronecker(d=3): Line 94  def kronecker(d=3):
94     @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
95     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
96     @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
97     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2.
98     """     """
99     return identityTensor(d)     return identityTensor(d)
100    
# Line 154  def identityTensor(d=3): Line 130  def identityTensor(d=3):
130     @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
131     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
132     @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
133     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
134     """     """
135     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
136         return escript.Data(identity((d.getDim(),)),d)         return escript.Data(identity((d.getDim(),)),d)
# Line 170  def identityTensor4(d=3): Line 146  def identityTensor4(d=3):
146     @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
147     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
148     @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
149     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
150     """     """
151     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
152         return escript.Data(identity((d.getDim(),d.getDim())),d)         return escript.Data(identity((d.getDim(),d.getDim())),d)
# Line 188  def unitVector(i=0,d=3): Line 164  def unitVector(i=0,d=3):
164     @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
165     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
166     @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
167     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 1     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
168     """     """
169     return kronecker(d)[i]     return kronecker(d)[i]
170    
# Line 244  def inf(arg): Line 220  def inf(arg):
220    
221      @param arg: argument      @param arg: argument
222      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
223      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
224      @rtype: C{float}      @rtype: C{float}
225      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
226      """      """
# Line 333  def commonDim(*args): Line 309  def commonDim(*args):
309      """      """
310      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.
311    
312      @param *args: given objects      @param args: given objects
313      @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
314               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
315      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 355  def testForZero(arg): Line 331  def testForZero(arg):
331    
332      @param arg: a given object      @param arg: a given object
333      @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}
334      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
335      @rtype : C{bool}      @rtype: C{bool}
336      """      """
337      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
338         return not Lsup(arg)>0.         return not Lsup(arg)>0.
# Line 459  def matchType(arg0=0.,arg1=0.): Line 435  def matchType(arg0=0.,arg1=0.):
435    
436  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
437      """      """
438            return representations of arg0 amd arg1 which ahve the same shape
   
     If shape is not given the shape "largest" shape of args is used.  
439    
440      @param args: a given ob      @param arg0: a given object
441      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
442      @return: True if the argument is identical to zero.      @param arg1: a given object
443      @rtype: C{list} of C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
444        @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
445        @rtype: C{tuple}
446      """      """
447      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
448      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
# Line 494  class Symbol(object): Line 470  class Symbol(object):
470         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
471         symbols or any other object.         symbols or any other object.
472    
473         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
474         @type arg: C{list}         @type args: C{list}
475         @param shape: the shape         @param shape: the shape
476         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
477         @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 538  class Symbol(object): Line 514  class Symbol(object):
514         """         """
515         the shape of the symbol.         the shape of the symbol.
516    
517         @return : the shape of the symbol.         @return: the shape of the symbol.
518         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
519         """         """
520         return self.__shape         return self.__shape
# Line 547  class Symbol(object): Line 523  class Symbol(object):
523         """         """
524         the spatial dimension         the spatial dimension
525    
526         @return : the spatial dimension         @return: the spatial dimension
527         @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.
528         """         """
529         return self.__dim         return self.__dim
# Line 571  class Symbol(object): Line 547  class Symbol(object):
547         """         """
548         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.
549    
550         @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].
551         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
552         @rtype: C{list} of objects         @rtype: C{list} of objects
553         @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.
554         """         """
555         out=[]         out=[]
556         for a in self.getArgument():         for a in self.getArgument():
# Line 698  class Symbol(object): Line 674  class Symbol(object):
674         """         """
675         returns -self.         returns -self.
676    
677         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
678         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
679         """         """
680         return self*(-1.)         return self*(-1.)
# Line 707  class Symbol(object): Line 683  class Symbol(object):
683         """         """
684         returns +self.         returns +self.
685    
686         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
687         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
688         """         """
689         return self*(1.)         return self*(1.)
690    
691     def __abs__(self):     def __abs__(self):
692         """         """
693         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
694         """         """
695         return Abs_Symbol(self)         return Abs_Symbol(self)
696    
# Line 724  class Symbol(object): Line 700  class Symbol(object):
700    
701         @param other: object to be added to this object         @param other: object to be added to this object
702         @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}.
703         @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}
704         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
705         """         """
706         return add(self,other)         return add(self,other)
# Line 735  class Symbol(object): Line 711  class Symbol(object):
711    
712         @param other: object this object is added to         @param other: object this object is added to
713         @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}.
714         @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
715         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
716         """         """
717         return add(other,self)         return add(other,self)
# Line 746  class Symbol(object): Line 722  class Symbol(object):
722    
723         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
724         @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}.
725         @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
726         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
727         """         """
728         return add(self,-other)         return add(self,-other)
# Line 757  class Symbol(object): Line 733  class Symbol(object):
733    
734         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
735         @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}.
736         @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}.
737         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
738         """         """
739         return add(-self,other)         return add(-self,other)
# Line 768  class Symbol(object): Line 744  class Symbol(object):
744    
745         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
746         @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}.
747         @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}.
748         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
749         """         """
750         return mult(self,other)         return mult(self,other)
# Line 779  class Symbol(object): Line 755  class Symbol(object):
755    
756         @param other: object this object is multiplied with         @param other: object this object is multiplied with
757         @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}.
758         @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.
759         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
760         """         """
761         return mult(other,self)         return mult(other,self)
# Line 790  class Symbol(object): Line 766  class Symbol(object):
766    
767         @param other: object dividing this object         @param other: object dividing this object
768         @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}.
769         @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}
770         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
771         """         """
772         return quotient(self,other)         return quotient(self,other)
# Line 801  class Symbol(object): Line 777  class Symbol(object):
777    
778         @param other: object dividing this object         @param other: object dividing this object
779         @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}.
780         @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
781         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
782         """         """
783         return quotient(other,self)         return quotient(other,self)
# Line 812  class Symbol(object): Line 788  class Symbol(object):
788    
789         @param other: exponent         @param other: exponent
790         @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}.
791         @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}
792         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
793         """         """
794         return power(self,other)         return power(self,other)
# Line 823  class Symbol(object): Line 799  class Symbol(object):
799    
800         @param other: basis         @param other: basis
801         @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}.
802         @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
803         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
804         """         """
805         return power(other,self)         return power(other,self)
# Line 834  class Symbol(object): Line 810  class Symbol(object):
810    
811         @param index: defines a         @param index: defines a
812         @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
813         @return: a S{Symbol} representing the slice defined by index         @return: a L{Symbol} representing the slice defined by index
814         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
815         """         """
816         return GetSlice_Symbol(self,index)         return GetSlice_Symbol(self,index)
# Line 844  class DependendSymbol(Symbol): Line 820  class DependendSymbol(Symbol):
820     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.
821     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  
822        
823     Example:     Example::
824        
825     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
826     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
827     print u1==u2       print u1==u2
828     False       False
829        
830        but       but::
831    
832     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
833     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
834     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
835     print u1==u2, u1==u3       print u1==u2, u1==u3
836     True False       True False
837    
838     @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.
839     """     """
# Line 947  class GetSlice_Symbol(DependendSymbol): Line 923  class GetSlice_Symbol(DependendSymbol):
923        @type format: C{str}        @type format: C{str}
924        @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.
925        @rtype: C{str}        @rtype: C{str}
926        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
927        """        """
928        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
929           return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])           return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
# Line 983  def log10(arg): Line 959  def log10(arg):
959    
960     @param arg: argument     @param arg: argument
961     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
962     @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.
963     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
964     """     """
965     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1005  def wherePositive(arg): Line 981  def wherePositive(arg):
981    
982     @param arg: argument     @param arg: argument
983     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
984     @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.
985     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
986     """     """
987     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1051  class WherePositive_Symbol(DependendSymb Line 1027  class WherePositive_Symbol(DependendSymb
1027        @type format: C{str}        @type format: C{str}
1028        @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.
1029        @rtype: C{str}        @rtype: C{str}
1030        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1031        """        """
1032        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1033            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1087  def whereNegative(arg): Line 1063  def whereNegative(arg):
1063    
1064     @param arg: argument     @param arg: argument
1065     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1066     @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.
1067     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1068     """     """
1069     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1133  class WhereNegative_Symbol(DependendSymb Line 1109  class WhereNegative_Symbol(DependendSymb
1109        @type format: C{str}        @type format: C{str}
1110        @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.
1111        @rtype: C{str}        @rtype: C{str}
1112        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1113        """        """
1114        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1115            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1169  def whereNonNegative(arg): Line 1145  def whereNonNegative(arg):
1145    
1146     @param arg: argument     @param arg: argument
1147     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1148     @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.
1149     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1150     """     """
1151     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1199  def whereNonPositive(arg): Line 1175  def whereNonPositive(arg):
1175    
1176     @param arg: argument     @param arg: argument
1177     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1178     @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.
1179     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1180     """     """
1181     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1231  def whereZero(arg,tol=0.): Line 1207  def whereZero(arg,tol=0.):
1207     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1208     @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.
1209     @type tol: C{float}     @type tol: C{float}
1210     @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.
1211     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1212     """     """
1213     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1239  def whereZero(arg,tol=0.): Line 1215  def whereZero(arg,tol=0.):
1215        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1216        return out        return out
1217     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1218        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1219     elif isinstance(arg,float):     elif isinstance(arg,float):
1220        if abs(arg)<=tol:        if abs(arg)<=tol:
1221          return 1.          return 1.
# Line 1280  class WhereZero_Symbol(DependendSymbol): Line 1253  class WhereZero_Symbol(DependendSymbol):
1253        @type format: C{str}        @type format: C{str}
1254        @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.
1255        @rtype: C{str}        @rtype: C{str}
1256        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1257        """        """
1258        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1259           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1314  def whereNonZero(arg,tol=0.): Line 1287  def whereNonZero(arg,tol=0.):
1287    
1288     @param arg: argument     @param arg: argument
1289     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1290     @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.
1291     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1292     """     """
1293     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1322  def whereNonZero(arg,tol=0.): Line 1295  def whereNonZero(arg,tol=0.):
1295        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1296        return out        return out
1297     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1298        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1299     elif isinstance(arg,float):     elif isinstance(arg,float):
1300        if abs(arg)>tol:        if abs(arg)>tol:
1301          return 1.          return 1.
# Line 1341  def whereNonZero(arg,tol=0.): Line 1311  def whereNonZero(arg,tol=0.):
1311     else:     else:
1312        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1313    
1314    def erf(arg):
1315       """
1316       returns erf of argument arg
1317    
1318       @param arg: argument
1319       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1320       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1321       @raises TypeError: if the type of the argument is not expected.
1322       """
1323       if isinstance(arg,escript.Data):
1324          return arg._erf()
1325       else:
1326          raise TypeError,"erf: Unknown argument type."
1327    
1328  def sin(arg):  def sin(arg):
1329     """     """
1330     returns sine of argument arg     returns sine of argument arg
1331    
1332     @param arg: argument     @param arg: argument
1333     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1334     @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.
1335     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1336     """     """
1337     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1385  class Sin_Symbol(DependendSymbol): Line 1369  class Sin_Symbol(DependendSymbol):
1369        @type format: C{str}        @type format: C{str}
1370        @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.
1371        @rtype: C{str}        @rtype: C{str}
1372        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1373        """        """
1374        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1375            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1437  def cos(arg): Line 1421  def cos(arg):
1421    
1422     @param arg: argument     @param arg: argument
1423     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1424     @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.
1425     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1426     """     """
1427     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1475  class Cos_Symbol(DependendSymbol): Line 1459  class Cos_Symbol(DependendSymbol):
1459        @type format: C{str}        @type format: C{str}
1460        @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.
1461        @rtype: C{str}        @rtype: C{str}
1462        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1463        """        """
1464        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1465            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1527  def tan(arg): Line 1511  def tan(arg):
1511    
1512     @param arg: argument     @param arg: argument
1513     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1514     @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.
1515     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1516     """     """
1517     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1565  class Tan_Symbol(DependendSymbol): Line 1549  class Tan_Symbol(DependendSymbol):
1549        @type format: C{str}        @type format: C{str}
1550        @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.
1551        @rtype: C{str}        @rtype: C{str}
1552        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1553        """        """
1554        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1555            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1617  def asin(arg): Line 1601  def asin(arg):
1601    
1602     @param arg: argument     @param arg: argument
1603     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1604     @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.
1605     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1606     """     """
1607     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1655  class Asin_Symbol(DependendSymbol): Line 1639  class Asin_Symbol(DependendSymbol):
1639        @type format: C{str}        @type format: C{str}
1640        @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.
1641        @rtype: C{str}        @rtype: C{str}
1642        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1643        """        """
1644        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1645            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1707  def acos(arg): Line 1691  def acos(arg):
1691    
1692     @param arg: argument     @param arg: argument
1693     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1694     @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.
1695     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1696     """     """
1697     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1745  class Acos_Symbol(DependendSymbol): Line 1729  class Acos_Symbol(DependendSymbol):
1729        @type format: C{str}        @type format: C{str}
1730        @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.
1731        @rtype: C{str}        @rtype: C{str}
1732        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1733        """        """
1734        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1735            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1797  def atan(arg): Line 1781  def atan(arg):
1781    
1782     @param arg: argument     @param arg: argument
1783     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1784     @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.
1785     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1786     """     """
1787     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1835  class Atan_Symbol(DependendSymbol): Line 1819  class Atan_Symbol(DependendSymbol):
1819        @type format: C{str}        @type format: C{str}
1820        @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.
1821        @rtype: C{str}        @rtype: C{str}
1822        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1823        """        """
1824        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1825            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1887  def sinh(arg): Line 1871  def sinh(arg):
1871    
1872     @param arg: argument     @param arg: argument
1873     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1874     @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.
1875     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1876     """     """
1877     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1925  class Sinh_Symbol(DependendSymbol): Line 1909  class Sinh_Symbol(DependendSymbol):
1909        @type format: C{str}        @type format: C{str}
1910        @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.
1911        @rtype: C{str}        @rtype: C{str}
1912        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1913        """        """
1914        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1915            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1977  def cosh(arg): Line 1961  def cosh(arg):
1961    
1962     @param arg: argument     @param arg: argument
1963     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1964     @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.
1965     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1966     """     """
1967     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2015  class Cosh_Symbol(DependendSymbol): Line 1999  class Cosh_Symbol(DependendSymbol):
1999        @type format: C{str}        @type format: C{str}
2000        @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.
2001        @rtype: C{str}        @rtype: C{str}
2002        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2003        """        """
2004        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2005            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2067  def tanh(arg): Line 2051  def tanh(arg):
2051    
2052     @param arg: argument     @param arg: argument
2053     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2054     @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.
2055     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2056     """     """
2057     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2105  class Tanh_Symbol(DependendSymbol): Line 2089  class Tanh_Symbol(DependendSymbol):
2089        @type format: C{str}        @type format: C{str}
2090        @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.
2091        @rtype: C{str}        @rtype: C{str}
2092        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2093        """        """
2094        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2095            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2157  def asinh(arg): Line 2141  def asinh(arg):
2141    
2142     @param arg: argument     @param arg: argument
2143     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2144     @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.
2145     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2146     """     """
2147     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2195  class Asinh_Symbol(DependendSymbol): Line 2179  class Asinh_Symbol(DependendSymbol):
2179        @type format: C{str}        @type format: C{str}
2180        @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.
2181        @rtype: C{str}        @rtype: C{str}
2182        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2183        """        """
2184        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2185            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2247  def acosh(arg): Line 2231  def acosh(arg):
2231    
2232     @param arg: argument     @param arg: argument
2233     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2234     @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.
2235     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2236     """     """
2237     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2285  class Acosh_Symbol(DependendSymbol): Line 2269  class Acosh_Symbol(DependendSymbol):
2269        @type format: C{str}        @type format: C{str}
2270        @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.
2271        @rtype: C{str}        @rtype: C{str}
2272        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2273        """        """
2274        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2275            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2337  def atanh(arg): Line 2321  def atanh(arg):
2321    
2322     @param arg: argument     @param arg: argument
2323     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2324     @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.
2325     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2326     """     """
2327     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2375  class Atanh_Symbol(DependendSymbol): Line 2359  class Atanh_Symbol(DependendSymbol):
2359        @type format: C{str}        @type format: C{str}
2360        @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.
2361        @rtype: C{str}        @rtype: C{str}
2362        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2363        """        """
2364        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2365            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2427  def exp(arg): Line 2411  def exp(arg):
2411    
2412     @param arg: argument     @param arg: argument
2413     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2414     @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.
2415     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2416     """     """
2417     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2465  class Exp_Symbol(DependendSymbol): Line 2449  class Exp_Symbol(DependendSymbol):
2449        @type format: C{str}        @type format: C{str}
2450        @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.
2451        @rtype: C{str}        @rtype: C{str}
2452        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2453        """        """
2454        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2455            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2517  def sqrt(arg): Line 2501  def sqrt(arg):
2501    
2502     @param arg: argument     @param arg: argument
2503     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2504     @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.
2505     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2506     """     """
2507     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2555  class Sqrt_Symbol(DependendSymbol): Line 2539  class Sqrt_Symbol(DependendSymbol):
2539        @type format: C{str}        @type format: C{str}
2540        @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.
2541        @rtype: C{str}        @rtype: C{str}
2542        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2543        """        """
2544        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2545            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2607  def log(arg): Line 2591  def log(arg):
2591    
2592     @param arg: argument     @param arg: argument
2593     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2594     @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.
2595     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2596     """     """
2597     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2645  class Log_Symbol(DependendSymbol): Line 2629  class Log_Symbol(DependendSymbol):
2629        @type format: C{str}        @type format: C{str}
2630        @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.
2631        @rtype: C{str}        @rtype: C{str}
2632        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2633        """        """
2634        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2635            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2697  def sign(arg): Line 2681  def sign(arg):
2681    
2682     @param arg: argument     @param arg: argument
2683     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2684     @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.
2685     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2686     """     """
2687     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2745  class Abs_Symbol(DependendSymbol): Line 2729  class Abs_Symbol(DependendSymbol):
2729        @type format: C{str}        @type format: C{str}
2730        @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.
2731        @rtype: C{str}        @rtype: C{str}
2732        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2733        """        """
2734        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2735            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2797  def minval(arg): Line 2781  def minval(arg):
2781    
2782     @param arg: argument     @param arg: argument
2783     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2784     @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.
2785     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2786     """     """
2787     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2838  class Minval_Symbol(DependendSymbol): Line 2822  class Minval_Symbol(DependendSymbol):
2822        @type format: C{str}        @type format: C{str}
2823        @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.
2824        @rtype: C{str}        @rtype: C{str}
2825        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2826        """        """
2827        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2828            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2874  def maxval(arg): Line 2858  def maxval(arg):
2858    
2859     @param arg: argument     @param arg: argument
2860     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2861     @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.
2862     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2863     """     """
2864     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2915  class Maxval_Symbol(DependendSymbol): Line 2899  class Maxval_Symbol(DependendSymbol):
2899        @type format: C{str}        @type format: C{str}
2900        @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.
2901        @rtype: C{str}        @rtype: C{str}
2902        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2903        """        """
2904        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2905            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2951  def length(arg): Line 2935  def length(arg):
2935    
2936     @param arg: argument     @param arg: argument
2937     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2938     @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.
2939     """     """
2940     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2941    
# Line 2961  def trace(arg,axis_offset=0): Line 2945  def trace(arg,axis_offset=0):
2945    
2946     @param arg: argument     @param arg: argument
2947     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2948     @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
2949                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
2950     @type axis_offset: C{int}     @type axis_offset: C{int}
2951     @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.
2952     @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 2970  def trace(arg,axis_offset=0): Line 2954  def trace(arg,axis_offset=0):
2954     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
2955        sh=arg.shape        sh=arg.shape
2956        if len(sh)<2:        if len(sh)<2:
2957          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2958        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
2959          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
2960        s1=1        s1=1
2961        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
2962        s2=1        s2=1
2963        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2964        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
2965          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)
2966        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))
2967        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
2968        for i1 in range(s1):        for i1 in range(s1):
# Line 2987  def trace(arg,axis_offset=0): Line 2971  def trace(arg,axis_offset=0):
2971        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2972        return out        return out
2973     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
2974        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
2975            raise ValueError,"rank of argument must be greater than 1"
2976          if axis_offset<0 or axis_offset>arg.getRank()-2:
2977            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
2978          s=list(arg.getShape())        
2979          if not s[axis_offset] == s[axis_offset+1]:
2980            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2981          return arg._trace(axis_offset)
2982     elif isinstance(arg,float):     elif isinstance(arg,float):
2983        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
2984     elif isinstance(arg,int):     elif isinstance(arg,int):
2985        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
2986     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
2987        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
2988     else:     else:
2989        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
2990    
 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  
2991  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
2992     """     """
2993     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 3045  class Trace_Symbol(DependendSymbol): Line 2997  class Trace_Symbol(DependendSymbol):
2997        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
2998        @param arg: argument of function        @param arg: argument of function
2999        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3000        @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
3001                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3002        @type axis_offset: C{int}        @type axis_offset: C{int}
3003        """        """
3004        if arg.getRank()<2:        if arg.getRank()<2:
3005          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3006        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3007          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
3008        s=list(arg.getShape())                s=list(arg.getShape())        
3009        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3010          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)
3011        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())
3012    
3013     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 3068  class Trace_Symbol(DependendSymbol): Line 3020  class Trace_Symbol(DependendSymbol):
3020        @type format: C{str}        @type format: C{str}
3021        @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.
3022        @rtype: C{str}        @rtype: C{str}
3023        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3024        """        """
3025        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3026           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3112  class Trace_Symbol(DependendSymbol): Line 3064  class Trace_Symbol(DependendSymbol):
3064    
3065  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3066     """     """
3067     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.
3068    
3069     @param arg: argument     @param arg: argument
3070     @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}
3071     @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.
3072                         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.
3073     @type axis_offset: C{int}     @type axis_offset: C{int}
3074     @return: transpose of arg     @return: transpose of arg
3075     @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 3126  def transpose(arg,axis_offset=None): Line 3078  def transpose(arg,axis_offset=None):
3078        if axis_offset==None: axis_offset=int(arg.rank/2)        if axis_offset==None: axis_offset=int(arg.rank/2)
3079        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))
3080     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3081        if axis_offset==None: axis_offset=int(arg.getRank()/2)        r=arg.getRank()
3082        return escript_transpose(arg,axis_offset)        if axis_offset==None: axis_offset=int(r/2)
3083          if axis_offset<0 or axis_offset>r:
3084            raise ValueError,"axis_offset must be between 0 and %s"%r
3085          return arg._transpose(axis_offset)
3086     elif isinstance(arg,float):     elif isinstance(arg,float):
3087        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3088          raise ValueError,"transpose: axis_offset must be 0 for float argument"          raise ValueError,"axis_offset must be 0 for float argument"
3089        return arg        return arg
3090     elif isinstance(arg,int):     elif isinstance(arg,int):
3091        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3092          raise ValueError,"transpose: axis_offset must be 0 for int argument"          raise ValueError,"axis_offset must be 0 for int argument"
3093        return float(arg)        return float(arg)
3094     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3095        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3096        return Transpose_Symbol(arg,axis_offset)        return Transpose_Symbol(arg,axis_offset)
3097     else:     else:
3098        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"Unknown argument type."
3099    
 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  
3100  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3101     """     """
3102     L{Symbol} representing the result of the transpose function     L{Symbol} representing the result of the transpose function
# Line 3216  class Transpose_Symbol(DependendSymbol): Line 3107  class Transpose_Symbol(DependendSymbol):
3107    
3108        @param arg: argument of function        @param arg: argument of function
3109        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3110         @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.
3111                         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.
3112        @type axis_offset: C{int}        @type axis_offset: C{int}
3113        """        """
3114        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3115        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3116          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()
3117        s=arg.getShape()        s=arg.getShape()
3118        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())
3119    
# Line 3236  class Transpose_Symbol(DependendSymbol): Line 3127  class Transpose_Symbol(DependendSymbol):
3127        @type format: C{str}        @type format: C{str}
3128        @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.
3129        @rtype: C{str}        @rtype: C{str}
3130        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3131        """        """
3132        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3133           return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3277  class Transpose_Symbol(DependendSymbol): Line 3168  class Transpose_Symbol(DependendSymbol):
3168           return identity(self.getShape())           return identity(self.getShape())
3169        else:        else:
3170           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3171    
3172    def swap_axes(arg,axis0=0,axis1=1):
3173       """
3174       returns the swap of arg by swaping the components axis0 and axis1
3175    
3176       @param arg: argument
3177       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3178       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3179       @type axis0: C{int}
3180       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3181       @type axis1: C{int}
3182       @return: C{arg} with swaped components
3183       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3184       """
3185       if axis0 > axis1:
3186          axis0,axis1=axis1,axis0
3187       if isinstance(arg,numarray.NumArray):
3188          return numarray.swapaxes(arg,axis0,axis1)
3189       elif isinstance(arg,escript.Data):
3190          return arg._swap_axes(axis0,axis1)
3191       elif isinstance(arg,float):
3192          raise TyepError,"float argument is not supported."
3193       elif isinstance(arg,int):
3194          raise TyepError,"int argument is not supported."
3195       elif isinstance(arg,Symbol):
3196          return SwapAxes_Symbol(arg,axis0,axis1)
3197       else:
3198          raise TypeError,"Unknown argument type."
3199    
3200    class SwapAxes_Symbol(DependendSymbol):
3201       """
3202       L{Symbol} representing the result of the swap function
3203       """
3204       def __init__(self,arg,axis0=0,axis1=1):
3205          """
3206          initialization of swap L{Symbol} with argument arg
3207    
3208          @param arg: argument
3209          @type arg: L{Symbol}.
3210          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3211          @type axis0: C{int}
3212          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3213          @type axis1: C{int}
3214          """
3215          if arg.getRank()<2:
3216             raise ValueError,"argument must have at least rank 2."
3217          if axis0<0 or axis0>arg.getRank()-1:
3218             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3219          if axis1<0 or axis1>arg.getRank()-1:
3220             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3221          if axis0 == axis1:
3222             raise ValueError,"axis indices must be different."
3223          if axis0 > axis1:
3224             axis0,axis1=axis1,axis0
3225          s=arg.getShape()
3226          s_out=[]
3227          for i in range(len(s)):
3228             if i == axis0:
3229                s_out.append(s[axis1])
3230             elif i == axis1:
3231                s_out.append(s[axis0])
3232             else:
3233                s_out.append(s[i])
3234          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3235    
3236       def getMyCode(self,argstrs,format="escript"):
3237          """
3238          returns a program code that can be used to evaluate the symbol.
3239    
3240          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3241          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3242          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3243          @type format: C{str}
3244          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3245          @rtype: C{str}
3246          @raise NotImplementedError: if the requested format is not available
3247          """
3248          if format=="escript" or format=="str"  or format=="text":
3249             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3250          else:
3251             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3252    
3253       def substitute(self,argvals):
3254          """
3255          assigns new values to symbols in the definition of the symbol.
3256          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3257    
3258          @param argvals: new values assigned to symbols
3259          @type argvals: C{dict} with keywords of type L{Symbol}.
3260          @return: result of the substitution process. Operations are executed as much as possible.
3261          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3262          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3263          """
3264          if argvals.has_key(self):
3265             arg=argvals[self]
3266             if self.isAppropriateValue(arg):
3267                return arg
3268             else:
3269                raise TypeError,"%s: new value is not appropriate."%str(self)
3270          else:
3271             arg=self.getSubstitutedArguments(argvals)
3272             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3273    
3274       def diff(self,arg):
3275          """
3276          differential of this object
3277    
3278          @param arg: the derivative is calculated with respect to arg
3279          @type arg: L{escript.Symbol}
3280          @return: derivative with respect to C{arg}
3281          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3282          """
3283          if arg==self:
3284             return identity(self.getShape())
3285          else:
3286             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3287    
3288  def symmetric(arg):  def symmetric(arg):
3289      """      """
3290      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 3289  def symmetric(arg): Line 3297  def symmetric(arg):
3297      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3298        if arg.rank==2:        if arg.rank==2:
3299          if not (arg.shape[0]==arg.shape[1]):          if not (arg.shape[0]==arg.shape[1]):
3300             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3301        elif arg.rank==4:        elif arg.rank==4:
3302          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]):
3303             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3304        else:        else:
3305          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3306        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3307      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3308        return escript_symmetric(arg)        if arg.getRank()==2:
3309            if not (arg.getShape()[0]==arg.getShape()[1]):
3310               raise ValueError,"argument must be square."
3311            return arg._symmetric()
3312          elif arg.getRank()==4:
3313            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3314               raise ValueError,"argument must be square."
3315            return arg._symmetric()
3316          else:
3317            raise ValueError,"rank 2 or 4 is required."
3318      elif isinstance(arg,float):      elif isinstance(arg,float):
3319        return arg        return arg
3320      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3305  def symmetric(arg): Line 3322  def symmetric(arg):
3322      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3323        if arg.getRank()==2:        if arg.getRank()==2:
3324          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3325             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3326        elif arg.getRank()==4:        elif arg.getRank()==4:
3327          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]):
3328             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3329        else:        else:
3330          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3331        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3332      else:      else:
3333        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3334    
 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  
   
3335  def nonsymmetric(arg):  def nonsymmetric(arg):
3336      """      """
3337      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 3356  def nonsymmetric(arg): Line 3352  def nonsymmetric(arg):
3352          raise ValueError,"nonsymmetric: rank 2 or 4 is required."          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3353        return (arg-transpose(arg))/2        return (arg-transpose(arg))/2
3354      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3355        return escript_nonsymmetric(arg)        if arg.getRank()==2:
3356            if not (arg.getShape()[0]==arg.getShape()[1]):
3357               raise ValueError,"argument must be square."
3358            return arg._nonsymmetric()
3359          elif arg.getRank()==4:
3360            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3361               raise ValueError,"argument must be square."
3362            return arg._nonsymmetric()
3363          else:
3364            raise ValueError,"rank 2 or 4 is required."
3365      elif isinstance(arg,float):      elif isinstance(arg,float):
3366        return arg        return arg
3367      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3374  def nonsymmetric(arg): Line 3379  def nonsymmetric(arg):
3379      else:      else:
3380        raise TypeError,"nonsymmetric: Unknown argument type."        raise TypeError,"nonsymmetric: Unknown argument type."
3381    
 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  
   
   
3382  def inverse(arg):  def inverse(arg):
3383      """      """
3384      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3385    
3386      @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.
3387      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3388      @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])
3389      @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
3390      @remark: for L{escript.Data} objects the dimension is restricted to 3.      @note: for L{escript.Data} objects the dimension is restricted to 3.
3391      """      """
3392        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3393      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3394        return numarray.linear_algebra.inverse(arg)        return numarray.linear_algebra.inverse(arg)
3395      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
# Line 3498  class Inverse_Symbol(DependendSymbol): Line 3482  class Inverse_Symbol(DependendSymbol):
3482        @type format: C{str}        @type format: C{str}
3483        @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.
3484        @rtype: C{str}        @rtype: C{str}
3485        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3486        """        """
3487        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3488           return "inverse(%s)"%argstrs[0]           return "inverse(%s)"%argstrs[0]
# Line 3538  class Inverse_Symbol(DependendSymbol): Line 3522  class Inverse_Symbol(DependendSymbol):
3522        if arg==self:        if arg==self:
3523           return identity(self.getShape())           return identity(self.getShape())
3524        else:        else:
3525           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3526    
3527  def eigenvalues(arg):  def eigenvalues(arg):
3528      """      """
# Line 3549  def eigenvalues(arg): Line 3533  def eigenvalues(arg):
3533      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3534      @return: the eigenvalues in increasing order.      @return: the eigenvalues in increasing order.
3535      @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.
3536      @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.
3537      """      """
3538      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3539        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
# Line 3619  def eigenvalues_and_eigenvectors(arg): Line 3603  def eigenvalues_and_eigenvectors(arg):
3603               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
3604               the eigenvector coresponding to the i-th eigenvalue.               the eigenvector coresponding to the i-th eigenvalue.
3605      @rtype: L{tuple} of L{escript.Data}.      @rtype: L{tuple} of L{escript.Data}.
3606      @remark: The dimension is restricted to 3.      @note: The dimension is restricted to 3.
3607      """      """
3608      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3609        raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"        raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
# Line 3692  class Add_Symbol(DependendSymbol): Line 3676  class Add_Symbol(DependendSymbol):
3676        @type format: C{str}        @type format: C{str}
3677        @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.
3678        @rtype: C{str}        @rtype: C{str}
3679        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3680        """        """
3681        if format=="str" or format=="text":        if format=="str" or format=="text":
3682           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 3791  class Mult_Symbol(DependendSymbol): Line 3775  class Mult_Symbol(DependendSymbol):
3775        @type format: C{str}        @type format: C{str}
3776        @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.
3777        @rtype: C{str}        @rtype: C{str}
3778        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3779        """        """
3780        if format=="str" or format=="text":        if format=="str" or format=="text":
3781           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3896  class Quotient_Symbol(DependendSymbol): Line 3880  class Quotient_Symbol(DependendSymbol):
3880        @type format: C{str}        @type format: C{str}
3881        @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.
3882        @rtype: C{str}        @rtype: C{str}
3883        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3884        """        """
3885        if format=="str" or format=="text":        if format=="str" or format=="text":
3886           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 4000  class Power_Symbol(DependendSymbol): Line 3984  class Power_Symbol(DependendSymbol):
3984        @type format: C{str}        @type format: C{str}
3985        @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.
3986        @rtype: C{str}        @rtype: C{str}
3987        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3988        """        """
3989        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
3990           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 4082  def minimum(*args): Line 4066  def minimum(*args):
4066            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4067      return out      return out
4068    
4069  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4070      """      """
4071      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4072    
4073      @param arg: argument      @param arg: argument
4074      @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}
4075      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4076      @type arg: C{float}      @type minval: C{float} or C{None}
4077      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4078      @type arg: C{float}      @type maxval: C{float} or C{None}
4079      @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
4080               less then maxval are unchanged.               less then maxval are unchanged.
4081      @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
4082      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4083      """      """
4084      if minval>maxval:      if not minval==None and not maxval==None:
4085         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4086      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4087        if minval == None:
4088            tmp=arg
4089        else:
4090            tmp=maximum(minval,arg)
4091        if maxval == None:
4092            return tmp
4093        else:
4094            return minimum(tmp,maxval)
4095    
4096        
4097  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 4116  def inner(arg0,arg1): Line 4108  def inner(arg0,arg1):
4108      @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}
4109      @param arg1: second argument      @param arg1: second argument
4110      @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}
4111      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4112      @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
4113      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4114      """      """
# Line 4126  def inner(arg0,arg1): Line 4118  def inner(arg0,arg1):
4118          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4119      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4120    
4121    def outer(arg0,arg1):
4122        """
4123        the outer product of the two argument:
4124    
4125        out[t,s]=arg0[t]*arg1[s]
4126    
4127        where
4128    
4129            - s runs through arg0.Shape
4130            - t runs through arg1.Shape
4131    
4132        @param arg0: first argument
4133        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4134        @param arg1: second argument
4135        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4136        @return: the outer product of arg0 and arg1 at each data point
4137        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4138        """
4139        return generalTensorProduct(arg0,arg1,axis_offset=0)
4140    
4141  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4142      """      """
4143        see L{matrix_mult}
4144        """
4145        return matrix_mult(arg0,arg1)
4146    
4147    def matrix_mult(arg0,arg1):
4148        """
4149      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4150    
4151      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4152    
4153            or      or
4154    
4155      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4156    
4157      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.
4158    
4159      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4160      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4154  def matrixmult(arg0,arg1): Line 4172  def matrixmult(arg0,arg1):
4172          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4173      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4174    
4175  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4176      """      """
4177      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  
4178      """      """
4179      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
4180    
4181    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4182      """      """
4183      the tensor product of the two argument:      the tensor product of the two argument:
   
4184            
4185      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4186    
4187      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4188    
4189                   or      or
4190    
4191      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4192    
# Line 4191  def tensormult(arg0,arg1): Line 4195  def tensormult(arg0,arg1):
4195    
4196      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]
4197                                
4198                   or      or
4199    
4200      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]
4201    
4202                   or      or
4203    
4204      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]
4205    
4206      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  
4207      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.
4208    
4209      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4210      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4216  def tensormult(arg0,arg1): Line 4220  def tensormult(arg0,arg1):
4220      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):
4221         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4222      else:      else:
4223          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4224    
4225  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4226      """      """
# Line 4224  def generalTensorProduct(arg0,arg1,axis_ Line 4228  def generalTensorProduct(arg0,arg1,axis_
4228    
4229      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4230    
4231      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:]  
4232    
4233      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]
4234      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4235            - t runs through arg1.Shape[axis_offset:]
4236    
4237      @param arg0: first argument      @param arg0: first argument
4238      @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}
4239      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4240      @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}
4241      @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.
4242      @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 4246  def generalTensorProduct(arg0,arg1,axis_ Line 4249  def generalTensorProduct(arg0,arg1,axis_
4249             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4250         else:         else:
4251             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4252                 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)
4253             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4254             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4255             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 4272  def generalTensorProduct(arg0,arg1,axis_ Line 4275  def generalTensorProduct(arg0,arg1,axis_
4275                                    
4276  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4277     """     """
4278     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4279     """     """
4280     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4281         """         """
4282         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4283    
4284         @param arg0: numerator         @param arg0: first argument
4285         @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}.
4286         @param arg1: denominator         @param arg1: second argument
4287         @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}.
4288         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4289         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4290         """         """
4291         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
# Line 4305  class GeneralTensorProduct_Symbol(Depend Line 4308  class GeneralTensorProduct_Symbol(Depend
4308        @type format: C{str}        @type format: C{str}
4309        @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.
4310        @rtype: C{str}        @rtype: C{str}
4311        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4312        """        """
4313        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4314           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 4333  class GeneralTensorProduct_Symbol(Depend Line 4336  class GeneralTensorProduct_Symbol(Depend
4336           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4337           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4338    
4339  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4340      "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!!!"
4341      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4342      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
4343      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  def transposed_matrix_mult(arg0,arg1):
4344      shape10=arg1.getShape()[:axis_offset]      """
4345      shape1=arg1.getShape()[axis_offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4346      if not shape01==shape10:  
4347          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]
4348    
4349      # whatr function space should be used? (this here is not good!)      or
4350      fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
4351      # create return value:      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4352      out=escript.Data(0.,tuple(shape0+shape1),fs)  
4353      #      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
     s0=[[]]  
     for k in shape0:  
           s=[]  
           for j in s0:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s0=s  
     s1=[[]]  
     for k in shape1:  
           s=[]  
           for j in s1:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s1=s  
     s01=[[]]  
     for k in shape01:  
           s=[]  
           for j in s01:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s01=s  
   
     for i0 in s0:  
        for i1 in s1:  
          s=escript.Scalar(0.,fs)  
          for i01 in s01:  
             s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))  
          out.__setitem__(tuple(i0+i1),s)  
     return out  
4354    
4355        The first dimension of arg0 and arg1 must match.
4356    
4357        @param arg0: first argument of rank 2
4358        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4359        @param arg1: second argument of at least rank 1
4360        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4361        @return: the product of the transposed of arg0 and arg1 at each data point
4362        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4363        @raise ValueError: if the shapes of the arguments are not appropriate
4364        """
4365        sh0=pokeShape(arg0)
4366        sh1=pokeShape(arg1)
4367        if not len(sh0)==2 :
4368            raise ValueError,"first argument must have rank 2"
4369        if not len(sh1)==2 and not len(sh1)==1:
4370            raise ValueError,"second argument must have rank 1 or 2"
4371        return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4372    
4373    def transposed_tensor_mult(arg0,arg1):
4374        """
4375        the tensor product of the transposed of the first and the second argument
4376        
4377        for arg0 of rank 2 this is
4378    
4379        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4380    
4381        or
4382    
4383        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4384    
4385      
4386        and for arg0 of rank 4 this is
4387    
4388        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4389                  
4390        or
4391    
4392        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4393    
4394        or
4395    
4396        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4397    
4398        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4399        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4400    
4401        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4402    
4403        @param arg0: first argument of rank 2 or 4
4404        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4405        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4406        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4407        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4408        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4409        """
4410        sh0=pokeShape(arg0)
4411        sh1=pokeShape(arg1)
4412        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4413           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4414        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4415           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4416        else:
4417            raise ValueError,"first argument must have rank 2 or 4"
4418    
4419    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4420        """
4421        generalized tensor product of transposed of arg0 and arg1:
4422    
4423        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4424    
4425        where
4426    
4427            - s runs through arg0.Shape[axis_offset:]
4428            - r runs trough arg0.Shape[:axis_offset]
4429            - t runs through arg1.Shape[axis_offset:]
4430    
4431        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4432        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4433    
4434        @param arg0: first argument
4435        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4436        @param arg1: second argument
4437        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4438        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4439        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4440        """
4441        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4442        arg0,arg1=matchType(arg0,arg1)
4443        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4444        if isinstance(arg0,numarray.NumArray):
4445           if isinstance(arg1,Symbol):
4446               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4447           else:
4448               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4449                   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)
4450               arg0_c=arg0.copy()
4451               arg1_c=arg1.copy()
4452               sh0,sh1=arg0.shape,arg1.shape
4453               d0,d1,d01=1,1,1
4454               for i in sh0[axis_offset:]: d0*=i
4455               for i in sh1[axis_offset:]: d1*=i
4456               for i in sh0[:axis_offset]: d01*=i
4457               arg0_c.resize((d01,d0))
4458               arg1_c.resize((d01,d1))
4459               out=numarray.zeros((d0,d1),numarray.Float64)
4460               for i0 in range(d0):
4461                        for i1 in range(d1):
4462                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4463               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4464               return out
4465        elif isinstance(arg0,escript.Data):
4466           if isinstance(arg1,Symbol):
4467               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4468           else:
4469               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4470        else:      
4471           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4472                    
4473    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4474       """
4475       Symbol representing the general tensor product of the transposed of arg0 and arg1
4476       """
4477       def __init__(self,arg0,arg1,axis_offset=0):
4478           """
4479           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4480    
4481           @param arg0: first argument
4482           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4483           @param arg1: second argument
4484           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4485           @raise ValueError: inconsistent dimensions of arguments.
4486           @note: if both arguments have a spatial dimension, they must equal.
4487           """
4488           sh_arg0=pokeShape(arg0)
4489           sh_arg1=pokeShape(arg1)
4490           sh01=sh_arg0[:axis_offset]
4491           sh10=sh_arg1[:axis_offset]
4492           sh0=sh_arg0[axis_offset:]
4493           sh1=sh_arg1[axis_offset:]
4494           if not sh01==sh10:
4495               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)
4496           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4497    
4498       def getMyCode(self,argstrs,format="escript"):
4499          """
4500          returns a program code that can be used to evaluate the symbol.
4501    
4502          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4503          @type argstrs: C{list} of length 2 of C{str}.
4504          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4505          @type format: C{str}
4506          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4507          @rtype: C{str}
4508          @raise NotImplementedError: if the requested format is not available
4509          """
4510          if format=="escript" or format=="str" or format=="text":
4511             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4512          else:
4513             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4514    
4515       def substitute(self,argvals):
4516          """
4517          assigns new values to symbols in the definition of the symbol.
4518          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4519    
4520          @param argvals: new values assigned to symbols
4521          @type argvals: C{dict} with keywords of type L{Symbol}.
4522          @return: result of the substitution process. Operations are executed as much as possible.
4523          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4524          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4525          """
4526          if argvals.has_key(self):
4527             arg=argvals[self]
4528             if self.isAppropriateValue(arg):
4529                return arg
4530             else:
4531                raise TypeError,"%s: new value is not appropriate."%str(self)
4532          else:
4533             args=self.getSubstitutedArguments(argvals)
4534             return generalTransposedTensorProduct(args[0],args[1],args[2])
4535    
4536    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4537        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4538        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4539    
4540    def matrix_transposed_mult(arg0,arg1):
4541        """
4542        matrix-transposed(matrix) product of the two argument:
4543    
4544        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4545    
4546        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4547    
4548        The last dimensions of arg0 and arg1 must match.
4549    
4550        @param arg0: first argument of rank 2
4551        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4552        @param arg1: second argument of rank 2
4553        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4554        @return: the product of arg0 and the transposed of arg1 at each data point
4555        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4556        @raise ValueError: if the shapes of the arguments are not appropriate
4557        """
4558        sh0=pokeShape(arg0)
4559        sh1=pokeShape(arg1)
4560        if not len(sh0)==2 :
4561            raise ValueError,"first argument must have rank 2"
4562        if not len(sh1)==2 and not len(sh1)==1:
4563            raise ValueError,"second argument must have rank 1 or 2"
4564        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4565    
4566    def tensor_transposed_mult(arg0,arg1):
4567        """
4568        the tensor product of the first and the transpose of the second argument
4569        
4570        for arg0 of rank 2 this is
4571    
4572        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4573    
4574        and for arg0 of rank 4 this is
4575    
4576        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4577                  
4578        or
4579    
4580        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4581    
4582        In the first case the the second dimension of arg0 and arg1 must match and  
4583        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4584    
4585        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4586    
4587        @param arg0: first argument of rank 2 or 4
4588        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4589        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4590        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4591        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4592        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4593        """
4594        sh0=pokeShape(arg0)
4595        sh1=pokeShape(arg1)
4596        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4597           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4598        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4599           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4600        else:
4601            raise ValueError,"first argument must have rank 2 or 4"
4602    
4603    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4604        """
4605        generalized tensor product of transposed of arg0 and arg1:
4606    
4607        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4608    
4609        where
4610    
4611            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4612            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4613            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4614    
4615        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4616        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4617    
4618        @param arg0: first argument
4619        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4620        @param arg1: second argument
4621        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4622        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4623        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4624        """
4625        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4626        arg0,arg1=matchType(arg0,arg1)
4627        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4628        if isinstance(arg0,numarray.NumArray):
4629           if isinstance(arg1,Symbol):
4630               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4631           else:
4632               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4633                   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)
4634               arg0_c=arg0.copy()
4635               arg1_c=arg1.copy()
4636               sh0,sh1=arg0.shape,arg1.shape
4637               d0,d1,d01=1,1,1
4638               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4639               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4640               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4641               arg0_c.resize((d0,d01))
4642               arg1_c.resize((d1,d01))
4643               out=numarray.zeros((d0,d1),numarray.Float64)
4644               for i0 in range(d0):
4645                        for i1 in range(d1):
4646                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4647               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4648               return out
4649        elif isinstance(arg0,escript.Data):
4650           if isinstance(arg1,Symbol):
4651               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4652           else:
4653               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4654        else:      
4655           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4656                    
4657    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4658       """
4659       Symbol representing the general tensor product of arg0 and the transpose of arg1
4660       """
4661       def __init__(self,arg0,arg1,axis_offset=0):
4662           """
4663           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4664    
4665           @param arg0: first argument
4666           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4667           @param arg1: second argument
4668           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4669           @raise ValueError: inconsistent dimensions of arguments.
4670           @note: if both arguments have a spatial dimension, they must equal.
4671           """
4672           sh_arg0=pokeShape(arg0)
4673           sh_arg1=pokeShape(arg1)
4674           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4675           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4676           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4677           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4678           if not sh01==sh10:
4679               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)
4680           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4681    
4682       def getMyCode(self,argstrs,format="escript"):
4683          """
4684          returns a program code that can be used to evaluate the symbol.
4685    
4686          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4687          @type argstrs: C{list} of length 2 of C{str}.
4688          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4689          @type format: C{str}
4690          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4691          @rtype: C{str}
4692          @raise NotImplementedError: if the requested format is not available
4693          """
4694          if format=="escript" or format=="str" or format=="text":
4695             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4696          else:
4697             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4698    
4699       def substitute(self,argvals):
4700          """
4701          assigns new values to symbols in the definition of the symbol.
4702          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4703    
4704          @param argvals: new values assigned to symbols
4705          @type argvals: C{dict} with keywords of type L{Symbol}.
4706          @return: result of the substitution process. Operations are executed as much as possible.
4707          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4708          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4709          """
4710          if argvals.has_key(self):
4711             arg=argvals[self]
4712             if self.isAppropriateValue(arg):
4713                return arg
4714             else:
4715                raise TypeError,"%s: new value is not appropriate."%str(self)
4716          else:
4717             args=self.getSubstitutedArguments(argvals)
4718             return generalTensorTransposedProduct(args[0],args[1],args[2])
4719    
4720    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4721        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4722        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4723    
4724  #=========================================================  #=========================================================
4725  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
# Line 4396  def grad(arg,where=None): Line 4741  def grad(arg,where=None):
4741                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4742      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4743      @return: gradient of arg.      @return: gradient of arg.
4744      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4745      """      """
4746      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4747         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 4436  class Grad_Symbol(DependendSymbol): Line 4781  class Grad_Symbol(DependendSymbol):
4781        @type format: C{str}        @type format: C{str}
4782        @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.
4783        @rtype: C{str}        @rtype: C{str}
4784        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4785        """        """
4786        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4787           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4489  def integrate(arg,where=None): Line 4834  def integrate(arg,where=None):
4834                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4835      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4836      @return: integral of arg.      @return: integral of arg.
4837      @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}      @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4838      """      """
4839      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4840         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
# Line 4527  class Integrate_Symbol(DependendSymbol): Line 4872  class Integrate_Symbol(DependendSymbol):
4872        @type format: C{str}        @type format: C{str}
4873        @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.
4874        @rtype: C{str}        @rtype: C{str}
4875        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4876        """        """
4877        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4878           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4579  def interpolate(arg,where): Line 4924  def interpolate(arg,where):
4924      @param where: FunctionSpace to be interpolated to      @param where: FunctionSpace to be interpolated to
4925      @type where: L{escript.FunctionSpace}      @type where: L{escript.FunctionSpace}
4926      @return: interpolated argument      @return: interpolated argument
4927      @rtype:  C{escript.Data} or L{Symbol}      @rtype: C{escript.Data} or L{Symbol}
4928      """      """
4929      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4930         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
# Line 4610  class Interpolate_Symbol(DependendSymbol Line 4955  class Interpolate_Symbol(DependendSymbol
4955        @type format: C{str}        @type format: C{str}
4956        @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.
4957        @rtype: C{str}        @rtype: C{str}
4958        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4959        """        """
4960        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4961           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4663  def div(arg,where=None): Line 5008  def div(arg,where=None):
5008                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
5009      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
5010      @return: divergence of arg.      @return: divergence of arg.
5011      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5012      """      """
5013      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5014          dim=arg.getDim()          dim=arg.getDim()
# Line 4685  def jump(arg,domain=None): Line 5030  def jump(arg,domain=None):
5030                     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.
5031      @type domain: C{None} or L{escript.Domain}      @type domain: C{None} or L{escript.Domain}
5032      @return: jump of arg      @return: jump of arg
5033      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5034      """      """
5035      if domain==None: domain=arg.getDomain()      if domain==None: domain=arg.getDomain()
5036      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
# Line 4697  def L2(arg): Line 5042  def L2(arg):
5042      @param arg: function which L2 to be calculated.      @param arg: function which L2 to be calculated.
5043      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
5044      @return: L2 norm of arg.      @return: L2 norm of arg.
5045      @rtype:  L{float} or L{Symbol}      @rtype: L{float} or L{Symbol}
5046      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5047      """      """
5048      return sqrt(integrate(inner(arg,arg)))      return sqrt(integrate(inner(arg,arg)))

Legend:
Removed from v.588  
changed lines
  Added in v.881

  ViewVC Help
Powered by ViewVC 1.1.26