/[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 587 by gross, Fri Mar 10 02:26:50 2006 UTC revision 785 by gross, Tue Jul 25 03:48:10 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    
 # 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  #   some helpers:  #   some helpers:
32  #=========================================================  #=========================================================
# Line 59  def saveVTK(filename,domain=None,**data) Line 34  def saveVTK(filename,domain=None,**data)
34      """      """
35      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.
36    
37      Example:      Example::
38    
39         tmp=Scalar(..)         tmp=Scalar(..)
40         v=Vector(..)         v=Vector(..)
# Line 87  def saveDX(filename,domain=None,**data): Line 62  def saveDX(filename,domain=None,**data):
62      """      """
63      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.
64    
65      Example:      Example::
66    
67         tmp=Scalar(..)         tmp=Scalar(..)
68         v=Vector(..)         v=Vector(..)
# Line 118  def kronecker(d=3): Line 93  def kronecker(d=3):
93     @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
94     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
95     @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
96     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2.
97     """     """
98     return identityTensor(d)     return identityTensor(d)
99    
# Line 154  def identityTensor(d=3): Line 129  def identityTensor(d=3):
129     @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
130     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
131     @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
132     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
133     """     """
134     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
135         return escript.Data(identity((d.getDim(),)),d)         return escript.Data(identity((d.getDim(),)),d)
# Line 170  def identityTensor4(d=3): Line 145  def identityTensor4(d=3):
145     @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
146     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
147     @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
148     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
149     """     """
150     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
151         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 163  def unitVector(i=0,d=3):
163     @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
164     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
165     @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
166     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 1     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
167     """     """
168     return kronecker(d)[i]     return kronecker(d)[i]
169    
# Line 244  def inf(arg): Line 219  def inf(arg):
219    
220      @param arg: argument      @param arg: argument
221      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
222      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
223      @rtype: C{float}      @rtype: C{float}
224      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
225      """      """
# Line 333  def commonDim(*args): Line 308  def commonDim(*args):
308      """      """
309      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.
310    
311      @param *args: given objects      @param args: given objects
312      @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
313               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
314      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 355  def testForZero(arg): Line 330  def testForZero(arg):
330    
331      @param arg: a given object      @param arg: a given object
332      @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}
333      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
334      @rtype : C{bool}      @rtype: C{bool}
335      """      """
336      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
337         return not Lsup(arg)>0.         return not Lsup(arg)>0.
# Line 459  def matchType(arg0=0.,arg1=0.): Line 434  def matchType(arg0=0.,arg1=0.):
434    
435  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
436      """      """
437            return representations of arg0 amd arg1 which ahve the same shape
   
     If shape is not given the shape "largest" shape of args is used.  
438    
439      @param args: a given ob      @param arg0: a given object
440      @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}
441      @return: True if the argument is identical to zero.      @param arg1: a given object
442      @rtype: C{list} of C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
443        @return: arg0 and arg1 where copies are returned when the shape has to be changed.
444        @rtype: C{tuple}
445      """      """
446      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
447      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
# Line 494  class Symbol(object): Line 469  class Symbol(object):
469         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
470         symbols or any other object.         symbols or any other object.
471    
472         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
473         @type arg: C{list}         @type args: C{list}
474         @param shape: the shape         @param shape: the shape
475         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
476         @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 513  class Symbol(object):
513         """         """
514         the shape of the symbol.         the shape of the symbol.
515    
516         @return : the shape of the symbol.         @return: the shape of the symbol.
517         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
518         """         """
519         return self.__shape         return self.__shape
# Line 547  class Symbol(object): Line 522  class Symbol(object):
522         """         """
523         the spatial dimension         the spatial dimension
524    
525         @return : the spatial dimension         @return: the spatial dimension
526         @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.
527         """         """
528         return self.__dim         return self.__dim
# Line 571  class Symbol(object): Line 546  class Symbol(object):
546         """         """
547         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.
548    
549         @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].
550         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
551         @rtype: C{list} of objects         @rtype: C{list} of objects
552         @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.
553         """         """
554         out=[]         out=[]
555         for a in self.getArgument():         for a in self.getArgument():
# Line 698  class Symbol(object): Line 673  class Symbol(object):
673         """         """
674         returns -self.         returns -self.
675    
676         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
677         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
678         """         """
679         return self*(-1.)         return self*(-1.)
# Line 707  class Symbol(object): Line 682  class Symbol(object):
682         """         """
683         returns +self.         returns +self.
684    
685         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
686         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
687         """         """
688         return self*(1.)         return self*(1.)
689    
690     def __abs__(self):     def __abs__(self):
691         """         """
692         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
693         """         """
694         return Abs_Symbol(self)         return Abs_Symbol(self)
695    
# Line 724  class Symbol(object): Line 699  class Symbol(object):
699    
700         @param other: object to be added to this object         @param other: object to be added to this object
701         @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}.
702         @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}
703         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
704         """         """
705         return add(self,other)         return add(self,other)
# Line 735  class Symbol(object): Line 710  class Symbol(object):
710    
711         @param other: object this object is added to         @param other: object this object is added to
712         @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}.
713         @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
714         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
715         """         """
716         return add(other,self)         return add(other,self)
# Line 746  class Symbol(object): Line 721  class Symbol(object):
721    
722         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
723         @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}.
724         @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
725         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
726         """         """
727         return add(self,-other)         return add(self,-other)
# Line 757  class Symbol(object): Line 732  class Symbol(object):
732    
733         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
734         @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}.
735         @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}.
736         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
737         """         """
738         return add(-self,other)         return add(-self,other)
# Line 768  class Symbol(object): Line 743  class Symbol(object):
743    
744         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
745         @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}.
746         @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}.
747         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
748         """         """
749         return mult(self,other)         return mult(self,other)
# Line 779  class Symbol(object): Line 754  class Symbol(object):
754    
755         @param other: object this object is multiplied with         @param other: object this object is multiplied with
756         @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}.
757         @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.
758         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
759         """         """
760         return mult(other,self)         return mult(other,self)
# Line 790  class Symbol(object): Line 765  class Symbol(object):
765    
766         @param other: object dividing this object         @param other: object dividing this object
767         @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}.
768         @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}
769         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
770         """         """
771         return quotient(self,other)         return quotient(self,other)
# Line 801  class Symbol(object): Line 776  class Symbol(object):
776    
777         @param other: object dividing this object         @param other: object dividing this object
778         @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}.
779         @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
780         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
781         """         """
782         return quotient(other,self)         return quotient(other,self)
# Line 812  class Symbol(object): Line 787  class Symbol(object):
787    
788         @param other: exponent         @param other: exponent
789         @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}.
790         @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}
791         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
792         """         """
793         return power(self,other)         return power(self,other)
# Line 823  class Symbol(object): Line 798  class Symbol(object):
798    
799         @param other: basis         @param other: basis
800         @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}.
801         @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
802         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
803         """         """
804         return power(other,self)         return power(other,self)
# Line 834  class Symbol(object): Line 809  class Symbol(object):
809    
810         @param index: defines a         @param index: defines a
811         @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
812         @return: a S{Symbol} representing the slice defined by index         @return: a L{Symbol} representing the slice defined by index
813         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
814         """         """
815         return GetSlice_Symbol(self,index)         return GetSlice_Symbol(self,index)
# Line 844  class DependendSymbol(Symbol): Line 819  class DependendSymbol(Symbol):
819     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.
820     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  
821        
822     Example:     Example::
823        
824     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
825     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
826     print u1==u2       print u1==u2
827     False       False
828        
829        but       but::
830    
831     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
832     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
833     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
834     print u1==u2, u1==u3       print u1==u2, u1==u3
835     True False       True False
836    
837     @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.
838     """     """
# Line 947  class GetSlice_Symbol(DependendSymbol): Line 922  class GetSlice_Symbol(DependendSymbol):
922        @type format: C{str}        @type format: C{str}
923        @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.
924        @rtype: C{str}        @rtype: C{str}
925        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
926        """        """
927        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
928           return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])           return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
# Line 983  def log10(arg): Line 958  def log10(arg):
958    
959     @param arg: argument     @param arg: argument
960     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
961     @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.
962     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
963     """     """
964     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1005  def wherePositive(arg): Line 980  def wherePositive(arg):
980    
981     @param arg: argument     @param arg: argument
982     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
983     @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.
984     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
985     """     """
986     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1051  class WherePositive_Symbol(DependendSymb Line 1026  class WherePositive_Symbol(DependendSymb
1026        @type format: C{str}        @type format: C{str}
1027        @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.
1028        @rtype: C{str}        @rtype: C{str}
1029        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1030        """        """
1031        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1032            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1087  def whereNegative(arg): Line 1062  def whereNegative(arg):
1062    
1063     @param arg: argument     @param arg: argument
1064     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1065     @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.
1066     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1067     """     """
1068     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1133  class WhereNegative_Symbol(DependendSymb Line 1108  class WhereNegative_Symbol(DependendSymb
1108        @type format: C{str}        @type format: C{str}
1109        @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.
1110        @rtype: C{str}        @rtype: C{str}
1111        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1112        """        """
1113        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1114            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1169  def whereNonNegative(arg): Line 1144  def whereNonNegative(arg):
1144    
1145     @param arg: argument     @param arg: argument
1146     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1147     @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.
1148     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1149     """     """
1150     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1199  def whereNonPositive(arg): Line 1174  def whereNonPositive(arg):
1174    
1175     @param arg: argument     @param arg: argument
1176     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1177     @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.
1178     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1179     """     """
1180     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1231  def whereZero(arg,tol=0.): Line 1206  def whereZero(arg,tol=0.):
1206     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1207     @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.
1208     @type tol: C{float}     @type tol: C{float}
1209     @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.
1210     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1211     """     """
1212     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1239  def whereZero(arg,tol=0.): Line 1214  def whereZero(arg,tol=0.):
1214        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1215        return out        return out
1216     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1217        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1218     elif isinstance(arg,float):     elif isinstance(arg,float):
1219        if abs(arg)<=tol:        if abs(arg)<=tol:
1220          return 1.          return 1.
# Line 1280  class WhereZero_Symbol(DependendSymbol): Line 1252  class WhereZero_Symbol(DependendSymbol):
1252        @type format: C{str}        @type format: C{str}
1253        @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.
1254        @rtype: C{str}        @rtype: C{str}
1255        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1256        """        """
1257        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1258           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 1286  def whereNonZero(arg,tol=0.):
1286    
1287     @param arg: argument     @param arg: argument
1288     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1289     @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.
1290     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1291     """     """
1292     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1322  def whereNonZero(arg,tol=0.): Line 1294  def whereNonZero(arg,tol=0.):
1294        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1295        return out        return out
1296     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1297        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1298     elif isinstance(arg,float):     elif isinstance(arg,float):
1299        if abs(arg)>tol:        if abs(arg)>tol:
1300          return 1.          return 1.
# Line 1347  def sin(arg): Line 1316  def sin(arg):
1316    
1317     @param arg: argument     @param arg: argument
1318     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1319     @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.
1320     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1321     """     """
1322     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1385  class Sin_Symbol(DependendSymbol): Line 1354  class Sin_Symbol(DependendSymbol):
1354        @type format: C{str}        @type format: C{str}
1355        @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.
1356        @rtype: C{str}        @rtype: C{str}
1357        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1358        """        """
1359        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1360            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1437  def cos(arg): Line 1406  def cos(arg):
1406    
1407     @param arg: argument     @param arg: argument
1408     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1409     @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.
1410     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1411     """     """
1412     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1475  class Cos_Symbol(DependendSymbol): Line 1444  class Cos_Symbol(DependendSymbol):
1444        @type format: C{str}        @type format: C{str}
1445        @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.
1446        @rtype: C{str}        @rtype: C{str}
1447        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1448        """        """
1449        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1450            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1527  def tan(arg): Line 1496  def tan(arg):
1496    
1497     @param arg: argument     @param arg: argument
1498     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1499     @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.
1500     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1501     """     """
1502     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1565  class Tan_Symbol(DependendSymbol): Line 1534  class Tan_Symbol(DependendSymbol):
1534        @type format: C{str}        @type format: C{str}
1535        @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.
1536        @rtype: C{str}        @rtype: C{str}
1537        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1538        """        """
1539        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1540            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1617  def asin(arg): Line 1586  def asin(arg):
1586    
1587     @param arg: argument     @param arg: argument
1588     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1589     @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.
1590     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1591     """     """
1592     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1655  class Asin_Symbol(DependendSymbol): Line 1624  class Asin_Symbol(DependendSymbol):
1624        @type format: C{str}        @type format: C{str}
1625        @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.
1626        @rtype: C{str}        @rtype: C{str}
1627        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1628        """        """
1629        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1630            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1707  def acos(arg): Line 1676  def acos(arg):
1676    
1677     @param arg: argument     @param arg: argument
1678     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1679     @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.
1680     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1681     """     """
1682     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1745  class Acos_Symbol(DependendSymbol): Line 1714  class Acos_Symbol(DependendSymbol):
1714        @type format: C{str}        @type format: C{str}
1715        @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.
1716        @rtype: C{str}        @rtype: C{str}
1717        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1718        """        """
1719        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1720            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1797  def atan(arg): Line 1766  def atan(arg):
1766    
1767     @param arg: argument     @param arg: argument
1768     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1769     @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.
1770     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1771     """     """
1772     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1835  class Atan_Symbol(DependendSymbol): Line 1804  class Atan_Symbol(DependendSymbol):
1804        @type format: C{str}        @type format: C{str}
1805        @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.
1806        @rtype: C{str}        @rtype: C{str}
1807        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1808        """        """
1809        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1810            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1887  def sinh(arg): Line 1856  def sinh(arg):
1856    
1857     @param arg: argument     @param arg: argument
1858     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1859     @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.
1860     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1861     """     """
1862     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1925  class Sinh_Symbol(DependendSymbol): Line 1894  class Sinh_Symbol(DependendSymbol):
1894        @type format: C{str}        @type format: C{str}
1895        @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.
1896        @rtype: C{str}        @rtype: C{str}
1897        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1898        """        """
1899        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1900            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1977  def cosh(arg): Line 1946  def cosh(arg):
1946    
1947     @param arg: argument     @param arg: argument
1948     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1949     @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.
1950     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1951     """     """
1952     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2015  class Cosh_Symbol(DependendSymbol): Line 1984  class Cosh_Symbol(DependendSymbol):
1984        @type format: C{str}        @type format: C{str}
1985        @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.
1986        @rtype: C{str}        @rtype: C{str}
1987        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1988        """        """
1989        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1990            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2067  def tanh(arg): Line 2036  def tanh(arg):
2036    
2037     @param arg: argument     @param arg: argument
2038     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2039     @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.
2040     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2041     """     """
2042     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2105  class Tanh_Symbol(DependendSymbol): Line 2074  class Tanh_Symbol(DependendSymbol):
2074        @type format: C{str}        @type format: C{str}
2075        @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.
2076        @rtype: C{str}        @rtype: C{str}
2077        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2078        """        """
2079        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2080            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2157  def asinh(arg): Line 2126  def asinh(arg):
2126    
2127     @param arg: argument     @param arg: argument
2128     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2129     @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.
2130     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2131     """     """
2132     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2195  class Asinh_Symbol(DependendSymbol): Line 2164  class Asinh_Symbol(DependendSymbol):
2164        @type format: C{str}        @type format: C{str}
2165        @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.
2166        @rtype: C{str}        @rtype: C{str}
2167        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2168        """        """
2169        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2170            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2247  def acosh(arg): Line 2216  def acosh(arg):
2216    
2217     @param arg: argument     @param arg: argument
2218     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2219     @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.
2220     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2221     """     """
2222     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2285  class Acosh_Symbol(DependendSymbol): Line 2254  class Acosh_Symbol(DependendSymbol):
2254        @type format: C{str}        @type format: C{str}
2255        @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.
2256        @rtype: C{str}        @rtype: C{str}
2257        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2258        """        """
2259        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2260            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2337  def atanh(arg): Line 2306  def atanh(arg):
2306    
2307     @param arg: argument     @param arg: argument
2308     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2309     @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.
2310     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2311     """     """
2312     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2375  class Atanh_Symbol(DependendSymbol): Line 2344  class Atanh_Symbol(DependendSymbol):
2344        @type format: C{str}        @type format: C{str}
2345        @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.
2346        @rtype: C{str}        @rtype: C{str}
2347        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2348        """        """
2349        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2350            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2427  def exp(arg): Line 2396  def exp(arg):
2396    
2397     @param arg: argument     @param arg: argument
2398     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2399     @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.
2400     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2401     """     """
2402     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2465  class Exp_Symbol(DependendSymbol): Line 2434  class Exp_Symbol(DependendSymbol):
2434        @type format: C{str}        @type format: C{str}
2435        @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.
2436        @rtype: C{str}        @rtype: C{str}
2437        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2438        """        """
2439        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2440            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2517  def sqrt(arg): Line 2486  def sqrt(arg):
2486    
2487     @param arg: argument     @param arg: argument
2488     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2489     @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.
2490     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2491     """     """
2492     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2555  class Sqrt_Symbol(DependendSymbol): Line 2524  class Sqrt_Symbol(DependendSymbol):
2524        @type format: C{str}        @type format: C{str}
2525        @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.
2526        @rtype: C{str}        @rtype: C{str}
2527        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2528        """        """
2529        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2530            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2607  def log(arg): Line 2576  def log(arg):
2576    
2577     @param arg: argument     @param arg: argument
2578     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2579     @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.
2580     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2581     """     """
2582     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2645  class Log_Symbol(DependendSymbol): Line 2614  class Log_Symbol(DependendSymbol):
2614        @type format: C{str}        @type format: C{str}
2615        @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.
2616        @rtype: C{str}        @rtype: C{str}
2617        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2618        """        """
2619        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2620            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2697  def sign(arg): Line 2666  def sign(arg):
2666    
2667     @param arg: argument     @param arg: argument
2668     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2669     @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.
2670     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2671     """     """
2672     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2745  class Abs_Symbol(DependendSymbol): Line 2714  class Abs_Symbol(DependendSymbol):
2714        @type format: C{str}        @type format: C{str}
2715        @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.
2716        @rtype: C{str}        @rtype: C{str}
2717        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2718        """        """
2719        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2720            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2797  def minval(arg): Line 2766  def minval(arg):
2766    
2767     @param arg: argument     @param arg: argument
2768     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2769     @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.
2770     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2771     """     """
2772     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2838  class Minval_Symbol(DependendSymbol): Line 2807  class Minval_Symbol(DependendSymbol):
2807        @type format: C{str}        @type format: C{str}
2808        @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.
2809        @rtype: C{str}        @rtype: C{str}
2810        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2811        """        """
2812        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2813            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2874  def maxval(arg): Line 2843  def maxval(arg):
2843    
2844     @param arg: argument     @param arg: argument
2845     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2846     @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.
2847     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2848     """     """
2849     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2915  class Maxval_Symbol(DependendSymbol): Line 2884  class Maxval_Symbol(DependendSymbol):
2884        @type format: C{str}        @type format: C{str}
2885        @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.
2886        @rtype: C{str}        @rtype: C{str}
2887        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2888        """        """
2889        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2890            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2951  def length(arg): Line 2920  def length(arg):
2920    
2921     @param arg: argument     @param arg: argument
2922     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2923     @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.
2924     """     """
2925     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2926    
# Line 2970  def trace(arg,axis_offset=0): Line 2939  def trace(arg,axis_offset=0):
2939     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
2940        sh=arg.shape        sh=arg.shape
2941        if len(sh)<2:        if len(sh)<2:
2942          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2943        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
2944          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
2945        s1=1        s1=1
2946        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
2947        s2=1        s2=1
2948        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2949        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
2950          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)
2951        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))
2952        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
2953        for i1 in range(s1):        for i1 in range(s1):
# Line 2987  def trace(arg,axis_offset=0): Line 2956  def trace(arg,axis_offset=0):
2956        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2957        return out        return out
2958     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
2959        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
2960            raise ValueError,"rank of argument must be greater than 1"
2961          if axis_offset<0 or axis_offset>arg.getRank()-2:
2962            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
2963          s=list(arg.getShape())        
2964          if not s[axis_offset] == s[axis_offset+1]:
2965            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2966          return arg._matrixtrace(axis_offset)
2967     elif isinstance(arg,float):     elif isinstance(arg,float):
2968        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
2969     elif isinstance(arg,int):     elif isinstance(arg,int):
2970        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
2971     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
2972        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
2973     else:     else:
2974        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
2975    
 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  
2976  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
2977     """     """
2978     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 3050  class Trace_Symbol(DependendSymbol): Line 2987  class Trace_Symbol(DependendSymbol):
2987        @type axis_offset: C{int}        @type axis_offset: C{int}
2988        """        """
2989        if arg.getRank()<2:        if arg.getRank()<2:
2990          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2991        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
2992          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
2993        s=list(arg.getShape())                s=list(arg.getShape())        
2994        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
2995          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)
2996        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())
2997    
2998     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 3068  class Trace_Symbol(DependendSymbol): Line 3005  class Trace_Symbol(DependendSymbol):
3005        @type format: C{str}        @type format: C{str}
3006        @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.
3007        @rtype: C{str}        @rtype: C{str}
3008        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3009        """        """
3010        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3011           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3126  def transpose(arg,axis_offset=None): Line 3063  def transpose(arg,axis_offset=None):
3063        if axis_offset==None: axis_offset=int(arg.rank/2)        if axis_offset==None: axis_offset=int(arg.rank/2)
3064        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))
3065     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3066        if axis_offset==None: axis_offset=int(arg.getRank()/2)        r=arg.getRank()
3067        return escript_transpose(arg,axis_offset)        if axis_offset==None: axis_offset=int(r/2)
3068          if axis_offset<0 or axis_offset>r:
3069            raise ValueError,"axis_offset must be between 0 and %s"%r
3070          return arg._transpose(axis_offset)
3071     elif isinstance(arg,float):     elif isinstance(arg,float):
3072        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3073          raise ValueError,"transpose: axis_offset must be 0 for float argument"          raise ValueError,"axis_offset must be 0 for float argument"
3074        return arg        return arg
3075     elif isinstance(arg,int):     elif isinstance(arg,int):
3076        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3077          raise ValueError,"transpose: axis_offset must be 0 for int argument"          raise ValueError,"axis_offset must be 0 for int argument"
3078        return float(arg)        return float(arg)
3079     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3080        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3081        return Transpose_Symbol(arg,axis_offset)        return Transpose_Symbol(arg,axis_offset)
3082     else:     else:
3083        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"Unknown argument type."
3084    
 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  
3085  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3086     """     """
3087     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 3092  class Transpose_Symbol(DependendSymbol):
3092    
3093        @param arg: argument of function        @param arg: argument of function
3094        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3095         @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 axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3096                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3097        @type axis_offset: C{int}        @type axis_offset: C{int}
3098        """        """
3099        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3100        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3101          raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%r
3102        s=arg.getShape()        s=arg.getShape()
3103        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())
3104    
# Line 3236  class Transpose_Symbol(DependendSymbol): Line 3112  class Transpose_Symbol(DependendSymbol):
3112        @type format: C{str}        @type format: C{str}
3113        @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.
3114        @rtype: C{str}        @rtype: C{str}
3115        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3116        """        """
3117        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3118           return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3289  def symmetric(arg): Line 3165  def symmetric(arg):
3165      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3166        if arg.rank==2:        if arg.rank==2:
3167          if not (arg.shape[0]==arg.shape[1]):          if not (arg.shape[0]==arg.shape[1]):
3168             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3169        elif arg.rank==4:        elif arg.rank==4:
3170          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]):
3171             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3172        else:        else:
3173          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3174        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3175      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3176        return escript_symmetric(arg)        if arg.getRank()==2:
3177            if not (arg.getShape()[0]==arg.getShape()[1]):
3178               raise ValueError,"argument must be square."
3179            return arg._symmetric()
3180          elif arg.getRank()==4:
3181            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3182               raise ValueError,"argument must be square."
3183            return arg._symmetric()
3184          else:
3185            raise ValueError,"rank 2 or 4 is required."
3186      elif isinstance(arg,float):      elif isinstance(arg,float):
3187        return arg        return arg
3188      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3305  def symmetric(arg): Line 3190  def symmetric(arg):
3190      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3191        if arg.getRank()==2:        if arg.getRank()==2:
3192          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3193             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3194        elif arg.getRank()==4:        elif arg.getRank()==4:
3195          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]):
3196             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3197        else:        else:
3198          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3199        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3200      else:      else:
3201        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3202    
 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  
   
3203  def nonsymmetric(arg):  def nonsymmetric(arg):
3204      """      """
3205      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 3220  def nonsymmetric(arg):
3220          raise ValueError,"nonsymmetric: rank 2 or 4 is required."          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3221        return (arg-transpose(arg))/2        return (arg-transpose(arg))/2
3222      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3223        return escript_nonsymmetric(arg)        if arg.getRank()==2:
3224            if not (arg.getShape()[0]==arg.getShape()[1]):
3225               raise ValueError,"argument must be square."
3226            return arg._nonsymmetric()
3227          elif arg.getRank()==4:
3228            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3229               raise ValueError,"argument must be square."
3230            return arg._nonsymmetric()
3231          else:
3232            raise ValueError,"rank 2 or 4 is required."
3233      elif isinstance(arg,float):      elif isinstance(arg,float):
3234        return arg        return arg
3235      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3374  def nonsymmetric(arg): Line 3247  def nonsymmetric(arg):
3247      else:      else:
3248        raise TypeError,"nonsymmetric: Unknown argument type."        raise TypeError,"nonsymmetric: Unknown argument type."
3249    
 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  
   
   
3250  def inverse(arg):  def inverse(arg):
3251      """      """
3252      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
# Line 3404  def inverse(arg): Line 3255  def inverse(arg):
3255      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3256      @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 matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3257      @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
3258      @remark: for L{escript.Data} objects the dimension is restricted to 3.      @note: for L{escript.Data} objects the dimension is restricted to 3.
3259      """      """
3260        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3261      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3262        return numarray.linear_algebra.inverse(arg)        return numarray.linear_algebra.inverse(arg)
3263      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
# Line 3498  class Inverse_Symbol(DependendSymbol): Line 3350  class Inverse_Symbol(DependendSymbol):
3350        @type format: C{str}        @type format: C{str}
3351        @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.
3352        @rtype: C{str}        @rtype: C{str}
3353        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3354        """        """
3355        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3356           return "inverse(%s)"%argstrs[0]           return "inverse(%s)"%argstrs[0]
# Line 3538  class Inverse_Symbol(DependendSymbol): Line 3390  class Inverse_Symbol(DependendSymbol):
3390        if arg==self:        if arg==self:
3391           return identity(self.getShape())           return identity(self.getShape())
3392        else:        else:
3393           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3394    
3395  def eigenvalues(arg):  def eigenvalues(arg):
3396      """      """
# Line 3549  def eigenvalues(arg): Line 3401  def eigenvalues(arg):
3401      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3402      @return: the eigenvalues in increasing order.      @return: the eigenvalues in increasing order.
3403      @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.
3404      @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.
3405      """      """
3406      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3407        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
# Line 3610  def eigenvalues(arg): Line 3462  def eigenvalues(arg):
3462    
3463  def eigenvalues_and_eigenvectors(arg):  def eigenvalues_and_eigenvectors(arg):
3464      """      """
3465      returns the eigenvalues of the square matrix arg.      returns the eigenvalues and eigenvectors of the square matrix arg.
3466    
3467      @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.
3468                  arg must be symmetric, ie. transpose(arg)==arg (this is not checked).                  arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3469      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{escript.Data}
3470      @return: the eigenvalues in increasing order.      @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3471      @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.               eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3472      @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.               the eigenvector coresponding to the i-th eigenvalue.
3473        @rtype: L{tuple} of L{escript.Data}.
3474        @note: The dimension is restricted to 3.
3475      """      """
3476      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3477        raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"        raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
# Line 3690  class Add_Symbol(DependendSymbol): Line 3544  class Add_Symbol(DependendSymbol):
3544        @type format: C{str}        @type format: C{str}
3545        @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.
3546        @rtype: C{str}        @rtype: C{str}
3547        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3548        """        """
3549        if format=="str" or format=="text":        if format=="str" or format=="text":
3550           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 3789  class Mult_Symbol(DependendSymbol): Line 3643  class Mult_Symbol(DependendSymbol):
3643        @type format: C{str}        @type format: C{str}
3644        @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.
3645        @rtype: C{str}        @rtype: C{str}
3646        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3647        """        """
3648        if format=="str" or format=="text":        if format=="str" or format=="text":
3649           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3894  class Quotient_Symbol(DependendSymbol): Line 3748  class Quotient_Symbol(DependendSymbol):
3748        @type format: C{str}        @type format: C{str}
3749        @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.
3750        @rtype: C{str}        @rtype: C{str}
3751        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3752        """        """
3753        if format=="str" or format=="text":        if format=="str" or format=="text":
3754           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3998  class Power_Symbol(DependendSymbol): Line 3852  class Power_Symbol(DependendSymbol):
3852        @type format: C{str}        @type format: C{str}
3853        @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.
3854        @rtype: C{str}        @rtype: C{str}
3855        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3856        """        """
3857        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
3858           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 4087  def clip(arg,minval=0.,maxval=1.): Line 3941  def clip(arg,minval=0.,maxval=1.):
3941      @param arg: argument      @param arg: argument
3942      @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}
3943      @param minval: lower range      @param minval: lower range
3944      @type arg: C{float}      @type minval: C{float}
3945      @param maxval: upper range      @param maxval: upper range
3946      @type arg: C{float}      @type maxval: C{float}
3947      @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
3948               less then maxval are unchanged.               less then maxval are unchanged.
3949      @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
# Line 4114  def inner(arg0,arg1): Line 3968  def inner(arg0,arg1):
3968      @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}
3969      @param arg1: second argument      @param arg1: second argument
3970      @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}
3971      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
3972      @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
3973      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
3974      """      """
# Line 4124  def inner(arg0,arg1): Line 3978  def inner(arg0,arg1):
3978          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
3979      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
3980    
3981  def matrixmult(arg0,arg1):  def outer(arg0,arg1):
3982        """
3983        the outer product of the two argument:
3984    
3985        out[t,s]=arg0[t]*arg1[s]
3986    
3987        where
3988    
3989            - s runs through arg0.Shape
3990            - t runs through arg1.Shape
3991    
3992        @param arg0: first argument
3993        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
3994        @param arg1: second argument
3995        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
3996        @return: the outer product of arg0 and arg1 at each data point
3997        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3998        """
3999        return generalTensorProduct(arg0,arg1,axis_offset=0)
4000    
4001    def matrixmul(arg0,arg1):
4002        """
4003        see L{matrix_mult}
4004        """
4005        return matrix_mult(arg0,arg1)
4006    
4007    def matrix_mult(arg0,arg1):
4008      """      """
4009      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4010    
4011      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4012    
4013            or      or
4014    
4015      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4016    
4017      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.
4018    
4019      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4020      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4152  def matrixmult(arg0,arg1): Line 4032  def matrixmult(arg0,arg1):
4032          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4033      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4034    
4035  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4036      """      """
4037      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  
4038      """      """
4039      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
4040    
4041    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4042      """      """
4043      the tensor product of the two argument:      the tensor product of the two argument:
   
4044            
4045      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4046    
4047      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4048    
4049                   or      or
4050    
4051      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4052    
# Line 4189  def tensormult(arg0,arg1): Line 4055  def tensormult(arg0,arg1):
4055    
4056      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]
4057                                
4058                   or      or
4059    
4060      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]
4061    
4062                   or      or
4063    
4064      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]
4065    
4066      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  
4067      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.
4068    
4069      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4070      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4214  def tensormult(arg0,arg1): Line 4080  def tensormult(arg0,arg1):
4080      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):
4081         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4082      else:      else:
4083          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4084    
4085  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4086      """      """
# Line 4222  def generalTensorProduct(arg0,arg1,axis_ Line 4088  def generalTensorProduct(arg0,arg1,axis_
4088    
4089      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4090    
4091      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:]  
4092    
4093      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]
4094      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4095            - t runs through arg1.Shape[axis_offset:]
4096    
4097      @param arg0: first argument      @param arg0: first argument
4098      @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}
4099      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4100      @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}
4101      @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.
4102      @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 4244  def generalTensorProduct(arg0,arg1,axis_ Line 4109  def generalTensorProduct(arg0,arg1,axis_
4109             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4110         else:         else:
4111             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4112                 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)
4113             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4114             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4115             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 4270  def generalTensorProduct(arg0,arg1,axis_ Line 4135  def generalTensorProduct(arg0,arg1,axis_
4135                                    
4136  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4137     """     """
4138     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4139     """     """
4140     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4141         """         """
4142         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4143    
4144         @param arg0: numerator         @param arg0: first argument
4145         @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}.
4146         @param arg1: denominator         @param arg1: second argument
4147         @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}.
4148         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4149         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4150         """         """
4151         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
# Line 4303  class GeneralTensorProduct_Symbol(Depend Line 4168  class GeneralTensorProduct_Symbol(Depend
4168        @type format: C{str}        @type format: C{str}
4169        @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.
4170        @rtype: C{str}        @rtype: C{str}
4171        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4172        """        """
4173        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4174           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 4374  def escript_generalTensorProduct(arg0,ar Line 4239  def escript_generalTensorProduct(arg0,ar
4239      return out      return out
4240    
4241    
4242    def transposed_matrix_mult(arg0,arg1):
4243        """
4244        transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4245    
4246        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
4247    
4248        or
4249    
4250        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4251    
4252        The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4253    
4254        The first dimension of arg0 and arg1 must match.
4255    
4256        @param arg0: first argument of rank 2
4257        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4258        @param arg1: second argument of at least rank 1
4259        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4260        @return: the product of the transposed of arg0 and arg1 at each data point
4261        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4262        @raise ValueError: if the shapes of the arguments are not appropriate
4263        """
4264        sh0=pokeShape(arg0)
4265        sh1=pokeShape(arg1)
4266        if not len(sh0)==2 :
4267            raise ValueError,"first argument must have rank 2"
4268        if not len(sh1)==2 and not len(sh1)==1:
4269            raise ValueError,"second argument must have rank 1 or 2"
4270        return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4271    
4272    def transposed_tensor_mult(arg0,arg1):
4273        """
4274        the tensor product of the transposed of the first and the second argument
4275        
4276        for arg0 of rank 2 this is
4277    
4278        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4279    
4280        or
4281    
4282        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4283    
4284      
4285        and for arg0 of rank 4 this is
4286    
4287        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4288                  
4289        or
4290    
4291        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4292    
4293        or
4294    
4295        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4296    
4297        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4298        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4299    
4300        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4301    
4302        @param arg0: first argument of rank 2 or 4
4303        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4304        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4305        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4306        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4307        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4308        """
4309        sh0=pokeShape(arg0)
4310        sh1=pokeShape(arg1)
4311        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4312           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4313        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4314           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4315        else:
4316            raise ValueError,"first argument must have rank 2 or 4"
4317    
4318    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4319        """
4320        generalized tensor product of transposed of arg0 and arg1:
4321    
4322        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4323    
4324        where
4325    
4326            - s runs through arg0.Shape[axis_offset:]
4327            - r runs trough arg0.Shape[:axis_offset]
4328            - t runs through arg1.Shape[axis_offset:]
4329    
4330        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4331        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4332    
4333        @param arg0: first argument
4334        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4335        @param arg1: second argument
4336        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4337        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4338        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4339        """
4340        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4341        arg0,arg1=matchType(arg0,arg1)
4342        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4343        if isinstance(arg0,numarray.NumArray):
4344           if isinstance(arg1,Symbol):
4345               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4346           else:
4347               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4348                   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)
4349               arg0_c=arg0.copy()
4350               arg1_c=arg1.copy()
4351               sh0,sh1=arg0.shape,arg1.shape
4352               d0,d1,d01=1,1,1
4353               for i in sh0[axis_offset:]: d0*=i
4354               for i in sh1[axis_offset:]: d1*=i
4355               for i in sh0[:axis_offset]: d01*=i
4356               arg0_c.resize((d01,d0))
4357               arg1_c.resize((d01,d1))
4358               out=numarray.zeros((d0,d1),numarray.Float64)
4359               for i0 in range(d0):
4360                        for i1 in range(d1):
4361                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4362               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4363               return out
4364        elif isinstance(arg0,escript.Data):
4365           if isinstance(arg1,Symbol):
4366               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4367           else:
4368               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4369        else:      
4370           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4371                    
4372    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4373       """
4374       Symbol representing the general tensor product of the transposed of arg0 and arg1
4375       """
4376       def __init__(self,arg0,arg1,axis_offset=0):
4377           """
4378           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4379    
4380           @param arg0: first argument
4381           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4382           @param arg1: second argument
4383           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4384           @raise ValueError: inconsistent dimensions of arguments.
4385           @note: if both arguments have a spatial dimension, they must equal.
4386           """
4387           sh_arg0=pokeShape(arg0)
4388           sh_arg1=pokeShape(arg1)
4389           sh01=sh_arg0[:axis_offset]
4390           sh10=sh_arg1[:axis_offset]
4391           sh0=sh_arg0[axis_offset:]
4392           sh1=sh_arg1[axis_offset:]
4393           if not sh01==sh10:
4394               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)
4395           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4396    
4397       def getMyCode(self,argstrs,format="escript"):
4398          """
4399          returns a program code that can be used to evaluate the symbol.
4400    
4401          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4402          @type argstrs: C{list} of length 2 of C{str}.
4403          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4404          @type format: C{str}
4405          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4406          @rtype: C{str}
4407          @raise NotImplementedError: if the requested format is not available
4408          """
4409          if format=="escript" or format=="str" or format=="text":
4410             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4411          else:
4412             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4413    
4414       def substitute(self,argvals):
4415          """
4416          assigns new values to symbols in the definition of the symbol.
4417          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4418    
4419          @param argvals: new values assigned to symbols
4420          @type argvals: C{dict} with keywords of type L{Symbol}.
4421          @return: result of the substitution process. Operations are executed as much as possible.
4422          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4423          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4424          """
4425          if argvals.has_key(self):
4426             arg=argvals[self]
4427             if self.isAppropriateValue(arg):
4428                return arg
4429             else:
4430                raise TypeError,"%s: new value is not appropriate."%str(self)
4431          else:
4432             args=self.getSubstitutedArguments(argvals)
4433             return generalTransposedTensorProduct(args[0],args[1],args[2])
4434    
4435    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4436        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4437        # calculate the return shape:
4438        shape01=arg0.getShape()[:axis_offset]
4439        shape10=arg1.getShape()[:axis_offset]
4440        shape0=arg0.getShape()[axis_offset:]
4441        shape1=arg1.getShape()[axis_offset:]
4442        if not shape01==shape10:
4443            raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4444    
4445        # whatr function space should be used? (this here is not good!)
4446        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4447        # create return value:
4448        out=escript.Data(0.,tuple(shape0+shape1),fs)
4449        #
4450        s0=[[]]
4451        for k in shape0:
4452              s=[]
4453              for j in s0:
4454                    for i in range(k): s.append(j+[slice(i,i)])
4455              s0=s
4456        s1=[[]]
4457        for k in shape1:
4458              s=[]
4459              for j in s1:
4460                    for i in range(k): s.append(j+[slice(i,i)])
4461              s1=s
4462        s01=[[]]
4463        for k in shape01:
4464              s=[]
4465              for j in s01:
4466                    for i in range(k): s.append(j+[slice(i,i)])
4467              s01=s
4468    
4469        for i0 in s0:
4470           for i1 in s1:
4471             s=escript.Scalar(0.,fs)
4472             for i01 in s01:
4473                s+=arg0.__getitem__(tuple(i01+i0))*arg1.__getitem__(tuple(i01+i1))
4474             out.__setitem__(tuple(i0+i1),s)
4475        return out
4476    
4477    
4478    def matrix_transposed_mult(arg0,arg1):
4479        """
4480        matrix-transposed(matrix) product of the two argument:
4481    
4482        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4483    
4484        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4485    
4486        The last dimensions of arg0 and arg1 must match.
4487    
4488        @param arg0: first argument of rank 2
4489        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4490        @param arg1: second argument of rank 2
4491        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4492        @return: the product of arg0 and the transposed of arg1 at each data point
4493        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4494        @raise ValueError: if the shapes of the arguments are not appropriate
4495        """
4496        sh0=pokeShape(arg0)
4497        sh1=pokeShape(arg1)
4498        if not len(sh0)==2 :
4499            raise ValueError,"first argument must have rank 2"
4500        if not len(sh1)==2 and not len(sh1)==1:
4501            raise ValueError,"second argument must have rank 1 or 2"
4502        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4503    
4504    def tensor_transposed_mult(arg0,arg1):
4505        """
4506        the tensor product of the first and the transpose of the second argument
4507        
4508        for arg0 of rank 2 this is
4509    
4510        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4511    
4512        and for arg0 of rank 4 this is
4513    
4514        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4515                  
4516        or
4517    
4518        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4519    
4520        In the first case the the second dimension of arg0 and arg1 must match and  
4521        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4522    
4523        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4524    
4525        @param arg0: first argument of rank 2 or 4
4526        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4527        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4528        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4529        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4530        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4531        """
4532        sh0=pokeShape(arg0)
4533        sh1=pokeShape(arg1)
4534        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4535           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4536        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4537           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4538        else:
4539            raise ValueError,"first argument must have rank 2 or 4"
4540    
4541    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4542        """
4543        generalized tensor product of transposed of arg0 and arg1:
4544    
4545        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4546    
4547        where
4548    
4549            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4550            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4551            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4552    
4553        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4554        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4555    
4556        @param arg0: first argument
4557        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4558        @param arg1: second argument
4559        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4560        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4561        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4562        """
4563        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4564        arg0,arg1=matchType(arg0,arg1)
4565        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4566        if isinstance(arg0,numarray.NumArray):
4567           if isinstance(arg1,Symbol):
4568               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4569           else:
4570               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4571                   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)
4572               arg0_c=arg0.copy()
4573               arg1_c=arg1.copy()
4574               sh0,sh1=arg0.shape,arg1.shape
4575               d0,d1,d01=1,1,1
4576               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4577               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4578               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4579               arg0_c.resize((d0,d01))
4580               arg1_c.resize((d1,d01))
4581               out=numarray.zeros((d0,d1),numarray.Float64)
4582               for i0 in range(d0):
4583                        for i1 in range(d1):
4584                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4585               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4586               return out
4587        elif isinstance(arg0,escript.Data):
4588           if isinstance(arg1,Symbol):
4589               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4590           else:
4591               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4592        else:      
4593           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4594                    
4595    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4596       """
4597       Symbol representing the general tensor product of arg0 and the transpose of arg1
4598       """
4599       def __init__(self,arg0,arg1,axis_offset=0):
4600           """
4601           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4602    
4603           @param arg0: first argument
4604           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4605           @param arg1: second argument
4606           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4607           @raise ValueError: inconsistent dimensions of arguments.
4608           @note: if both arguments have a spatial dimension, they must equal.
4609           """
4610           sh_arg0=pokeShape(arg0)
4611           sh_arg1=pokeShape(arg1)
4612           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4613           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4614           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4615           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4616           if not sh01==sh10:
4617               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)
4618           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4619    
4620       def getMyCode(self,argstrs,format="escript"):
4621          """
4622          returns a program code that can be used to evaluate the symbol.
4623    
4624          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4625          @type argstrs: C{list} of length 2 of C{str}.
4626          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4627          @type format: C{str}
4628          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4629          @rtype: C{str}
4630          @raise NotImplementedError: if the requested format is not available
4631          """
4632          if format=="escript" or format=="str" or format=="text":
4633             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4634          else:
4635             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4636    
4637       def substitute(self,argvals):
4638          """
4639          assigns new values to symbols in the definition of the symbol.
4640          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4641    
4642          @param argvals: new values assigned to symbols
4643          @type argvals: C{dict} with keywords of type L{Symbol}.
4644          @return: result of the substitution process. Operations are executed as much as possible.
4645          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4646          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4647          """
4648          if argvals.has_key(self):
4649             arg=argvals[self]
4650             if self.isAppropriateValue(arg):
4651                return arg
4652             else:
4653                raise TypeError,"%s: new value is not appropriate."%str(self)
4654          else:
4655             args=self.getSubstitutedArguments(argvals)
4656             return generalTensorTransposedProduct(args[0],args[1],args[2])
4657    
4658    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4659        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4660        # calculate the return shape:
4661        shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4662        shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4663        shape10=arg1.getShape()[arg1.getRank()-axis_offset:]
4664        shape1=arg1.getShape()[:arg1.getRank()-axis_offset]
4665        if not shape01==shape10:
4666            raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4667    
4668        # whatr function space should be used? (this here is not good!)
4669        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4670        # create return value:
4671        out=escript.Data(0.,tuple(shape0+shape1),fs)
4672        #
4673        s0=[[]]
4674        for k in shape0:
4675              s=[]
4676              for j in s0:
4677                    for i in range(k): s.append(j+[slice(i,i)])
4678              s0=s
4679        s1=[[]]
4680        for k in shape1:
4681              s=[]
4682              for j in s1:
4683                    for i in range(k): s.append(j+[slice(i,i)])
4684              s1=s
4685        s01=[[]]
4686        for k in shape01:
4687              s=[]
4688              for j in s01:
4689                    for i in range(k): s.append(j+[slice(i,i)])
4690              s01=s
4691    
4692        for i0 in s0:
4693           for i1 in s1:
4694             s=escript.Scalar(0.,fs)
4695             for i01 in s01:
4696                s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i1+i01))
4697             out.__setitem__(tuple(i0+i1),s)
4698        return out
4699    
4700    
4701  #=========================================================  #=========================================================
4702  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
4703  #=========================================================  #=========================================================
# Line 4394  def grad(arg,where=None): Line 4718  def grad(arg,where=None):
4718                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4719      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4720      @return: gradient of arg.      @return: gradient of arg.
4721      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4722      """      """
4723      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4724         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 4434  class Grad_Symbol(DependendSymbol): Line 4758  class Grad_Symbol(DependendSymbol):
4758        @type format: C{str}        @type format: C{str}
4759        @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.
4760        @rtype: C{str}        @rtype: C{str}
4761        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4762        """        """
4763        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4764           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4487  def integrate(arg,where=None): Line 4811  def integrate(arg,where=None):
4811                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4812      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4813      @return: integral of arg.      @return: integral of arg.
4814      @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}      @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4815      """      """
4816      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4817         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
# Line 4525  class Integrate_Symbol(DependendSymbol): Line 4849  class Integrate_Symbol(DependendSymbol):
4849        @type format: C{str}        @type format: C{str}
4850        @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.
4851        @rtype: C{str}        @rtype: C{str}
4852        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4853        """        """
4854        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4855           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4577  def interpolate(arg,where): Line 4901  def interpolate(arg,where):
4901      @param where: FunctionSpace to be interpolated to      @param where: FunctionSpace to be interpolated to
4902      @type where: L{escript.FunctionSpace}      @type where: L{escript.FunctionSpace}
4903      @return: interpolated argument      @return: interpolated argument
4904      @rtype:  C{escript.Data} or L{Symbol}      @rtype: C{escript.Data} or L{Symbol}
4905      """      """
4906      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4907         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
# Line 4608  class Interpolate_Symbol(DependendSymbol Line 4932  class Interpolate_Symbol(DependendSymbol
4932        @type format: C{str}        @type format: C{str}
4933        @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.
4934        @rtype: C{str}        @rtype: C{str}
4935        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4936        """        """
4937        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4938           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4661  def div(arg,where=None): Line 4985  def div(arg,where=None):
4985                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4986      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4987      @return: divergence of arg.      @return: divergence of arg.
4988      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4989      """      """
4990      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4991          dim=arg.getDim()          dim=arg.getDim()
# Line 4683  def jump(arg,domain=None): Line 5007  def jump(arg,domain=None):
5007                     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.
5008      @type domain: C{None} or L{escript.Domain}      @type domain: C{None} or L{escript.Domain}
5009      @return: jump of arg      @return: jump of arg
5010      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5011      """      """
5012      if domain==None: domain=arg.getDomain()      if domain==None: domain=arg.getDomain()
5013      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
# Line 4695  def L2(arg): Line 5019  def L2(arg):
5019      @param arg: function which L2 to be calculated.      @param arg: function which L2 to be calculated.
5020      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
5021      @return: L2 norm of arg.      @return: L2 norm of arg.
5022      @rtype:  L{float} or L{Symbol}      @rtype: L{float} or L{Symbol}
5023      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5024      """      """
5025      return sqrt(integrate(inner(arg,arg)))      return sqrt(integrate(inner(arg,arg)))

Legend:
Removed from v.587  
changed lines
  Added in v.785

  ViewVC Help
Powered by ViewVC 1.1.26