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

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

Parent Directory | Revision Log | View Patch Patch

-revision 587 by gross,
Fri Mar 10 02:26:50 2006 UTC
+revision 804 by gross,
Thu Aug 10 01:12:16 2006 UTC
 Line 1
  # $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.
- #
  """
  Utility functions for escript
- @remark:  This module is under construction and is still tested!!!
  @var __author__: name of author
- @var __licence__: licence agreement
+ @var __copyright__: copyrights
+ @var __license__: licence agreement
  @var __url__: url entry point on documentation
  @var __version__: version
  @var __date__: date of the version
  """
  __author__="Lutz Gross, l.gross@uq.edu.au"
- __licence__="contact: esys@access.uq.edu.au"
+ __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
+                     http://www.access.edu.au
+                 Primary Business: Queensland, Australia"""
+ __license__="""Licensed under the Open Software License version 3.0
+              http://www.opensource.org/licenses/osl-3.0.php"""
  __url__="http://www.iservo.edu.au/esys/escript"
  __version__="$Revision$"
  __date__="$Date$"
-Line 30 
 __date__="$Date$"
+Line 24 
 __date__="$Date$"
  import math
  import numarray
- import numarray.linear_algebra
  import escript
  import os
- # 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
  #=========================================================
  #   some helpers:
  #=========================================================
-Line 59 
 def saveVTK(filename,domain=None,**data)
+Line 34 
 def saveVTK(filename,domain=None,**data)
      """
      writes a L{Data} objects into a files using the the VTK XML file format.
-     Example:
+     Example::
         tmp=Scalar(..)
         v=Vector(..)
-Line 87 
 def saveDX(filename,domain=None,**data):
+Line 62 
 def saveDX(filename,domain=None,**data):
      """
      writes a L{Data} objects into a files using the the DX file format.
-     Example:
+     Example::
         tmp=Scalar(..)
         v=Vector(..)
-Line 118 
 def kronecker(d=3):
+Line 93 
 def kronecker(d=3):
     @param d: dimension or an object that has the C{getDim} method defining the dimension
     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
-    @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2.
+    @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2.
     """
     return identityTensor(d)
-Line 154 
 def identityTensor(d=3):
+Line 129 
 def identityTensor(d=3):
     @param d: dimension or an object that has the C{getDim} method defining the dimension
     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
-    @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2
+    @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
     """
     if isinstance(d,escript.FunctionSpace):
         return escript.Data(identity((d.getDim(),)),d)
-Line 170 
 def identityTensor4(d=3):
+Line 145 
 def identityTensor4(d=3):
     @param d: dimension or an object that has the C{getDim} method defining the dimension
     @type d: C{int} or any object with a C{getDim} method
     @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
-    @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 4.
+    @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
     """
     if isinstance(d,escript.FunctionSpace):
         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):
     @param d: dimension or an object that has the C{getDim} method defining the dimension
     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
     @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise
-    @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 1
+    @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
     """
     return kronecker(d)[i]
-Line 244 
 def inf(arg):
+Line 219 
 def inf(arg):
      @param arg: argument
      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
-     @return : minimum value of arg over all components and all data points
+     @return: minimum value of arg over all components and all data points
      @rtype: C{float}
      @raise TypeError: if type of arg cannot be processed
      """
-Line 333 
 def commonDim(*args):
+Line 308 
 def commonDim(*args):
      """
      identifies, if possible, the spatial dimension across a set of objects which may or my not have a spatial dimension.
-     @param *args: given objects
+     @param args: given objects
      @return: the spatial dimension of the objects with identifiable dimension (see L{pokeDim}). If none the objects has
               a spatial dimension C{None} is returned.
      @rtype: C{int} or C{None}
-Line 355 
 def testForZero(arg):
+Line 330 
 def testForZero(arg):
      @param arg: a given object
      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}
-     @return : True if the argument is identical to zero.
+     @return: True if the argument is identical to zero.
-     @rtype : C{bool}
+     @rtype: C{bool}
      """
      if isinstance(arg,numarray.NumArray):
         return not Lsup(arg)>0.
-Line 459 
 def matchType(arg0=0.,arg1=0.):
+Line 434 
 def matchType(arg0=0.,arg1=0.):
  def matchShape(arg0,arg1):
      """
+     return representations of arg0 amd arg1 which ahve the same shape
-     If shape is not given the shape "largest" shape of args is used.
+     @param arg0: a given object
+     @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
-     @param args: a given ob
+     @param arg1: a given object
-     @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}
+     @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
-     @return: True if the argument is identical to zero.
+     @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
-     @rtype: C{list} of C{int}
+     @rtype: C{tuple}
      """
      sh=commonShape(arg0,arg1)
      sh0=pokeShape(arg0)
-Line 494 
 class Symbol(object):
+Line 469 
 class Symbol(object):
         Creates an instance of a symbol of a given shape. The symbol may depending on a list of arguments args which may be
         symbols or any other object.
-        @param arg: the arguments of the symbol.
+        @param args: the arguments of the symbol.
-        @type arg: C{list}
+        @type args: C{list}
         @param shape: the shape
         @type shape: C{tuple} of C{int}
         @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):
         """
         the shape of the symbol.
-        @return : the shape of the symbol.
+        @return: the shape of the symbol.
         @rtype: C{tuple} of C{int}
         """
         return self.__shape
-Line 547 
 class Symbol(object):
+Line 522 
 class Symbol(object):
         """
         the spatial dimension
-        @return : the spatial dimension
+        @return: the spatial dimension
         @rtype: C{int} if the dimension is defined. Otherwise C{None} is returned.
         """
         return self.__dim
-Line 571 
 class Symbol(object):
+Line 546 
 class Symbol(object):
         """
         substitutes symbols in the arguments of this object and returns the result as a list.
-        @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].
         @type argvals: C{dict} with keywords of type L{Symbol}.
         @rtype: C{list} of objects
-        @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.
         """
         out=[]
         for a in self.getArgument():
-Line 698 
 class Symbol(object):
+Line 673 
 class Symbol(object):
         """
         returns -self.
-        @return:  a S{Symbol} representing the negative of the object
+        @return:  a L{Symbol} representing the negative of the object
         @rtype: L{DependendSymbol}
         """
         return self*(-1.)
-Line 707 
 class Symbol(object):
+Line 682 
 class Symbol(object):
         """
         returns +self.
-        @return:  a S{Symbol} representing the positive of the object
+        @return:  a L{Symbol} representing the positive of the object
         @rtype: L{DependendSymbol}
         """
         return self*(1.)
     def __abs__(self):
         """
-        returns a S{Symbol} representing the absolute value of the object.
+        returns a L{Symbol} representing the absolute value of the object.
         """
         return Abs_Symbol(self)
-Line 724 
 class Symbol(object):
+Line 699 
 class Symbol(object):
         @param other: object to be added to this object
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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}
         @rtype: L{DependendSymbol}
         """
         return add(self,other)
-Line 735 
 class Symbol(object):
+Line 710 
 class Symbol(object):
         @param other: object this object is added to
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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
         @rtype: L{DependendSymbol}
         """
         return add(other,self)
-Line 746 
 class Symbol(object):
+Line 721 
 class Symbol(object):
         @param other: object to be subtracted from this object
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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
         @rtype: L{DependendSymbol}
         """
         return add(self,-other)
-Line 757 
 class Symbol(object):
+Line 732 
 class Symbol(object):
         @param other: object this object is been subtracted from
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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}.
         @rtype: L{DependendSymbol}
         """
         return add(-self,other)
-Line 768 
 class Symbol(object):
+Line 743 
 class Symbol(object):
         @param other: object to be mutiplied by this object
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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}.
         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
         """
         return mult(self,other)
-Line 779 
 class Symbol(object):
+Line 754 
 class Symbol(object):
         @param other: object this object is multiplied with
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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.
         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
         """
         return mult(other,self)
-Line 790 
 class Symbol(object):
+Line 765 
 class Symbol(object):
         @param other: object dividing this object
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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}
         @rtype: L{DependendSymbol}
         """
         return quotient(self,other)
-Line 801 
 class Symbol(object):
+Line 776 
 class Symbol(object):
         @param other: object dividing this object
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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
         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
         """
         return quotient(other,self)
-Line 812 
 class Symbol(object):
+Line 787 
 class Symbol(object):
         @param other: exponent
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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}
         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
         """
         return power(self,other)
-Line 823 
 class Symbol(object):
+Line 798 
 class Symbol(object):
         @param other: basis
         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @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
         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
         """
         return power(other,self)
-Line 834 
 class Symbol(object):
+Line 809 
 class Symbol(object):
         @param index: defines a
         @type index: C{slice} or C{int} or a C{tuple} of them
-        @return: a S{Symbol} representing the slice defined by index
+        @return: a L{Symbol} representing the slice defined by index
         @rtype: L{DependendSymbol}
         """
         return GetSlice_Symbol(self,index)
-Line 844 
 class DependendSymbol(Symbol):
+Line 819 
 class DependendSymbol(Symbol):
     DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.
     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
-    Example:
+    Example::
-    u1=Symbol(shape=(3,4),dim=2,args=[4.])
+      u1=Symbol(shape=(3,4),dim=2,args=[4.])
-    u2=Symbol(shape=(3,4),dim=2,args=[4.])
+      u2=Symbol(shape=(3,4),dim=2,args=[4.])
-    print u1==u2
+      print u1==u2
-    False
+      False
-       but
+    but::
-    u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
+      u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
-    u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
+      u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
-    u3=DependendSymbol(shape=(2,),dim=2,args=[4.])
+      u3=DependendSymbol(shape=(2,),dim=2,args=[4.])
-    print u1==u2, u1==u3
+      print u1==u2, u1==u3
-    True False
+      True False
     @note: DependendSymbol should be used as return value of functions with L{Symbol} arguments. This will allow the optimizer to remove redundant function calls.
     """
-Line 947 
 class GetSlice_Symbol(DependendSymbol):
+Line 922 
 class GetSlice_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str"  or format=="text":
           return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
-Line 983 
 def log10(arg):
+Line 958 
 def log10(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1005 
 def wherePositive(arg):
+Line 980 
 def wherePositive(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1051 
 class WherePositive_Symbol(DependendSymb
+Line 1026 
 class WherePositive_Symbol(DependendSymb
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 1087 
 def whereNegative(arg):
+Line 1062 
 def whereNegative(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1133 
 class WhereNegative_Symbol(DependendSymb
+Line 1108 
 class WhereNegative_Symbol(DependendSymb
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 1169 
 def whereNonNegative(arg):
+Line 1144 
 def whereNonNegative(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1199 
 def whereNonPositive(arg):
+Line 1174 
 def whereNonPositive(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1231 
 def whereZero(arg,tol=0.):
+Line 1206 
 def whereZero(arg,tol=0.):
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
     @param tol: tolerance. values with absolute value less then tol are accepted as zero.
     @type tol: C{float}
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1239 
 def whereZero(arg,tol=0.):
+Line 1214 
 def whereZero(arg,tol=0.):
        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
        return out
     elif isinstance(arg,escript.Data):
-       if tol>0.:
+       return arg._whereZero(tol)
-          return whereNegative(abs(arg)-tol)
-       else:
-          return arg._whereZero()
     elif isinstance(arg,float):
        if abs(arg)<=tol:
          return 1.
-Line 1280 
 class WhereZero_Symbol(DependendSymbol):
+Line 1252 
 class WhereZero_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str"  or format=="text":
           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
-Line 1314 
 def whereNonZero(arg,tol=0.):
+Line 1286 
 def whereNonZero(arg,tol=0.):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1322 
 def whereNonZero(arg,tol=0.):
+Line 1294 
 def whereNonZero(arg,tol=0.):
        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
        return out
     elif isinstance(arg,escript.Data):
-       if tol>0.:
+       return arg._whereNonZero(tol)
-          return 1.-whereZero(arg,tol)
-       else:
-          return arg._whereNonZero()
     elif isinstance(arg,float):
        if abs(arg)>tol:
          return 1.
-Line 1347 
 def sin(arg):
+Line 1316 
 def sin(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1385 
 class Sin_Symbol(DependendSymbol):
+Line 1354 
 class Sin_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 1437 
 def cos(arg):
+Line 1406 
 def cos(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1475 
 class Cos_Symbol(DependendSymbol):
+Line 1444 
 class Cos_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 1527 
 def tan(arg):
+Line 1496 
 def tan(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1565 
 class Tan_Symbol(DependendSymbol):
+Line 1534 
 class Tan_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 1617 
 def asin(arg):
+Line 1586 
 def asin(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1655 
 class Asin_Symbol(DependendSymbol):
+Line 1624 
 class Asin_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 1707 
 def acos(arg):
+Line 1676 
 def acos(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1745 
 class Acos_Symbol(DependendSymbol):
+Line 1714 
 class Acos_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 1797 
 def atan(arg):
+Line 1766 
 def atan(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1835 
 class Atan_Symbol(DependendSymbol):
+Line 1804 
 class Atan_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 1887 
 def sinh(arg):
+Line 1856 
 def sinh(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 1925 
 class Sinh_Symbol(DependendSymbol):
+Line 1894 
 class Sinh_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 1977 
 def cosh(arg):
+Line 1946 
 def cosh(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2015 
 class Cosh_Symbol(DependendSymbol):
+Line 1984 
 class Cosh_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2067 
 def tanh(arg):
+Line 2036 
 def tanh(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2105 
 class Tanh_Symbol(DependendSymbol):
+Line 2074 
 class Tanh_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2157 
 def asinh(arg):
+Line 2126 
 def asinh(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2195 
 class Asinh_Symbol(DependendSymbol):
+Line 2164 
 class Asinh_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2247 
 def acosh(arg):
+Line 2216 
 def acosh(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2285 
 class Acosh_Symbol(DependendSymbol):
+Line 2254 
 class Acosh_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2337 
 def atanh(arg):
+Line 2306 
 def atanh(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2375 
 class Atanh_Symbol(DependendSymbol):
+Line 2344 
 class Atanh_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2427 
 def exp(arg):
+Line 2396 
 def exp(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2465 
 class Exp_Symbol(DependendSymbol):
+Line 2434 
 class Exp_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2517 
 def sqrt(arg):
+Line 2486 
 def sqrt(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2555 
 class Sqrt_Symbol(DependendSymbol):
+Line 2524 
 class Sqrt_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2607 
 def log(arg):
+Line 2576 
 def log(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2645 
 class Log_Symbol(DependendSymbol):
+Line 2614 
 class Log_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2697 
 def sign(arg):
+Line 2666 
 def sign(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2745 
 class Abs_Symbol(DependendSymbol):
+Line 2714 
 class Abs_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2797 
 def minval(arg):
+Line 2766 
 def minval(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2838 
 class Minval_Symbol(DependendSymbol):
+Line 2807 
 class Minval_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2874 
 def maxval(arg):
+Line 2843 
 def maxval(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     @raises TypeError: if the type of the argument is not expected.
     """
     if isinstance(arg,numarray.NumArray):
-Line 2915 
 class Maxval_Symbol(DependendSymbol):
+Line 2884 
 class Maxval_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if isinstance(argstrs,list):
            argstrs=argstrs[0]
-Line 2951 
 def length(arg):
+Line 2920 
 def length(arg):
     @param arg: argument
     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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.
     """
     return sqrt(inner(arg,arg))
-Line 2961 
 def trace(arg,axis_offset=0):
+Line 2930 
 def trace(arg,axis_offset=0):
     @param arg: argument
     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
+    @param axis_offset: C{axis_offset} to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
-                   axis_offset and axis_offset+1 must be equal.
+                   C{axis_offset} and axis_offset+1 must be equal.
     @type axis_offset: C{int}
     @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
-Line 2970 
 def trace(arg,axis_offset=0):
+Line 2939 
 def trace(arg,axis_offset=0):
     if isinstance(arg,numarray.NumArray):
        sh=arg.shape
        if len(sh)<2:
-         raise ValueError,"trace: rank of argument must be greater than 1"
+         raise ValueError,"rank of argument must be greater than 1"
        if axis_offset<0 or axis_offset>len(sh)-2:
-         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
        s1=1
        for i in range(axis_offset): s1*=sh[i]
        s2=1
        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
        if not sh[axis_offset] == sh[axis_offset+1]:
-         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)
        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
        out=numarray.zeros([s1,s2],numarray.Float64)
        for i1 in range(s1):
-Line 2987 
 def trace(arg,axis_offset=0):
+Line 2956 
 def trace(arg,axis_offset=0):
        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
        return out
     elif isinstance(arg,escript.Data):
-       return escript_trace(arg,axis_offset)
+       if arg.getRank()<2:
+         raise ValueError,"rank of argument must be greater than 1"
+       if axis_offset<0 or axis_offset>arg.getRank()-2:
+         raise ValueError,"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,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
+       return arg._trace(axis_offset)
     elif isinstance(arg,float):
-       raise TypeError,"trace: illegal argument type float."
+       raise TypeError,"illegal argument type float."
     elif isinstance(arg,int):
-       raise TypeError,"trace: illegal argument type int."
+       raise TypeError,"illegal argument type int."
     elif isinstance(arg,Symbol):
        return Trace_Symbol(arg,axis_offset)
     else:
-       raise TypeError,"trace: Unknown argument type."
+       raise TypeError,"Unknown argument type."
- 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
  class Trace_Symbol(DependendSymbol):
     """
     L{Symbol} representing the result of the trace function
-Line 3045 
 class Trace_Symbol(DependendSymbol):
+Line 2982 
 class Trace_Symbol(DependendSymbol):
        initialization of trace L{Symbol} with argument arg
        @param arg: argument of function
        @type arg: L{Symbol}.
-       @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
+       @param axis_offset: C{axis_offset} to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
-                   axis_offset and axis_offset+1 must be equal.
+                   C{axis_offset} and axis_offset+1 must be equal.
        @type axis_offset: C{int}
        """
        if arg.getRank()<2:
-         raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"
+         raise ValueError,"rank of argument must be greater than 1"
        if axis_offset<0 or axis_offset>arg.getRank()-2:
-         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
        s=list(arg.getShape())
        if not s[axis_offset] == s[axis_offset+1]:
-         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)
        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
     def getMyCode(self,argstrs,format="escript"):
-Line 3068 
 class Trace_Symbol(DependendSymbol):
+Line 3005 
 class Trace_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str"  or format=="text":
           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
-Line 3112 
 class Trace_Symbol(DependendSymbol):
+Line 3049 
 class Trace_Symbol(DependendSymbol):
  def transpose(arg,axis_offset=None):
     """
-    returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.
+    returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
     @param arg: argument
     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
-    @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
+    @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
-                        if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
+                        if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
     @type axis_offset: C{int}
     @return: transpose of arg
     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
-Line 3126 
 def transpose(arg,axis_offset=None):
+Line 3063 
 def transpose(arg,axis_offset=None):
        if axis_offset==None: axis_offset=int(arg.rank/2)
        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
     elif isinstance(arg,escript.Data):
-       if axis_offset==None: axis_offset=int(arg.getRank()/2)
+       r=arg.getRank()
-       return escript_transpose(arg,axis_offset)
+       if axis_offset==None: axis_offset=int(r/2)
+       if axis_offset<0 or axis_offset>r:
+         raise ValueError,"axis_offset must be between 0 and %s"%r
+       return arg._transpose(axis_offset)
     elif isinstance(arg,float):
        if not ( axis_offset==0 or axis_offset==None):
-         raise ValueError,"transpose: axis_offset must be 0 for float argument"
+         raise ValueError,"axis_offset must be 0 for float argument"
        return arg
     elif isinstance(arg,int):
        if not ( axis_offset==0 or axis_offset==None):
-         raise ValueError,"transpose: axis_offset must be 0 for int argument"
+         raise ValueError,"axis_offset must be 0 for int argument"
        return float(arg)
     elif isinstance(arg,Symbol):
        if axis_offset==None: axis_offset=int(arg.getRank()/2)
        return Transpose_Symbol(arg,axis_offset)
     else:
-       raise TypeError,"transpose: Unknown argument type."
+       raise TypeError,"Unknown argument type."
- 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
  class Transpose_Symbol(DependendSymbol):
     """
     L{Symbol} representing the result of the transpose function
-Line 3216 
 class Transpose_Symbol(DependendSymbol):
+Line 3092 
 class Transpose_Symbol(DependendSymbol):
        @param arg: argument of function
        @type arg: L{Symbol}.
-        @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
+       @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
-                        if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
+                        if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
        @type axis_offset: C{int}
        """
        if axis_offset==None: axis_offset=int(arg.getRank()/2)
        if axis_offset<0 or axis_offset>arg.getRank():
-         raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
+         raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
        s=arg.getShape()
        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
-Line 3236 
 class Transpose_Symbol(DependendSymbol):
+Line 3112 
 class Transpose_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str"  or format=="text":
           return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
-Line 3277 
 class Transpose_Symbol(DependendSymbol):
+Line 3153 
 class Transpose_Symbol(DependendSymbol):
           return identity(self.getShape())
        else:
           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
+ def swap_axes(arg,axis0=0,axis1=1):
+    """
+    returns the swap of arg by swaping the components axis0 and axis1
+    @param arg: argument
+    @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
+    @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
+    @type axis0: C{int}
+    @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
+    @type axis1: C{int}
+    @return: C{arg} with swaped components
+    @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
+    """
+    if axis0 > axis1:
+       axis0,axis1=axis1,axis0
+    if isinstance(arg,numarray.NumArray):
+       return numarray.swapaxes(arg,axis0,axis1)
+    elif isinstance(arg,escript.Data):
+       return arg._swap_axes(axis0,axis1)
+    elif isinstance(arg,float):
+       raise TyepError,"float argument is not supported."
+    elif isinstance(arg,int):
+       raise TyepError,"int argument is not supported."
+    elif isinstance(arg,Symbol):
+       return SwapAxes_Symbol(arg,axis0,axis1)
+    else:
+       raise TypeError,"Unknown argument type."
+ class SwapAxes_Symbol(DependendSymbol):
+    """
+    L{Symbol} representing the result of the swap function
+    """
+    def __init__(self,arg,axis0=0,axis1=1):
+       """
+       initialization of swap L{Symbol} with argument arg
+       @param arg: argument
+       @type arg: L{Symbol}.
+       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
+       @type axis0: C{int}
+       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
+       @type axis1: C{int}
+       """
+       if arg.getRank()<2:
+          raise ValueError,"argument must have at least rank 2."
+       if axis0<0 or axis0>arg.getRank()-1:
+          raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
+       if axis1<0 or axis1>arg.getRank()-1:
+          raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
+       if axis0 == axis1:
+          raise ValueError,"axis indices must be different."
+       if axis0 > axis1:
+          axis0,axis1=axis1,axis0
+       s=arg.getShape()
+       s_out=[]
+       for i in range(len(s)):
+          if i == axis0:
+             s_out.append(s[axis1])
+          elif i == axis1:
+             s_out.append(s[axis0])
+          else:
+             s_out.append(s[i])
+       super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
+    def getMyCode(self,argstrs,format="escript"):
+       """
+       returns a program code that can be used to evaluate the symbol.
+       @param argstrs: gives for each argument a string representing the argument for the evaluation.
+       @type argstrs: C{str} or a C{list} of length 1 of C{str}.
+       @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
+       @type format: C{str}
+       @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
+       @rtype: C{str}
+       @raise NotImplementedError: if the requested format is not available
+       """
+       if format=="escript" or format=="str"  or format=="text":
+          return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
+       else:
+          raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
+    def substitute(self,argvals):
+       """
+       assigns new values to symbols in the definition of the symbol.
+       The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
+       @param argvals: new values assigned to symbols
+       @type argvals: C{dict} with keywords of type L{Symbol}.
+       @return: result of the substitution process. Operations are executed as much as possible.
+       @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
+       @raise TypeError: if a value for a L{Symbol} cannot be substituted.
+       """
+       if argvals.has_key(self):
+          arg=argvals[self]
+          if self.isAppropriateValue(arg):
+             return arg
+          else:
+             raise TypeError,"%s: new value is not appropriate."%str(self)
+       else:
+          arg=self.getSubstitutedArguments(argvals)
+          return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
+    def diff(self,arg):
+       """
+       differential of this object
+       @param arg: the derivative is calculated with respect to arg
+       @type arg: L{escript.Symbol}
+       @return: derivative with respect to C{arg}
+       @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
+       """
+       if arg==self:
+          return identity(self.getShape())
+       else:
+          return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
  def symmetric(arg):
      """
      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
-Line 3289 
 def symmetric(arg):
+Line 3282 
 def symmetric(arg):
      if isinstance(arg,numarray.NumArray):
        if arg.rank==2:
          if not (arg.shape[0]==arg.shape[1]):
-            raise ValueError,"symmetric: argument must be square."
+            raise ValueError,"argument must be square."
        elif arg.rank==4:
          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
-            raise ValueError,"symmetric: argument must be square."
+            raise ValueError,"argument must be square."
        else:
-         raise ValueError,"symmetric: rank 2 or 4 is required."
+         raise ValueError,"rank 2 or 4 is required."
        return (arg+transpose(arg))/2
      elif isinstance(arg,escript.Data):
-       return escript_symmetric(arg)
+       if arg.getRank()==2:
+         if not (arg.getShape()[0]==arg.getShape()[1]):
+            raise ValueError,"argument must be square."
+         return arg._symmetric()
+       elif arg.getRank()==4:
+         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
+            raise ValueError,"argument must be square."
+         return arg._symmetric()
+       else:
+         raise ValueError,"rank 2 or 4 is required."
      elif isinstance(arg,float):
        return arg
      elif isinstance(arg,int):
-Line 3305 
 def symmetric(arg):
+Line 3307 
 def symmetric(arg):
      elif isinstance(arg,Symbol):
        if arg.getRank()==2:
          if not (arg.getShape()[0]==arg.getShape()[1]):
-            raise ValueError,"symmetric: argument must be square."
+            raise ValueError,"argument must be square."
        elif arg.getRank()==4:
          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
-            raise ValueError,"symmetric: argument must be square."
+            raise ValueError,"argument must be square."
        else:
-         raise ValueError,"symmetric: rank 2 or 4 is required."
+         raise ValueError,"rank 2 or 4 is required."
        return (arg+transpose(arg))/2
      else:
        raise TypeError,"symmetric: Unknown argument type."
- 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
  def nonsymmetric(arg):
      """
      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
-Line 3356 
 def nonsymmetric(arg):
+Line 3337 
 def nonsymmetric(arg):
          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
        return (arg-transpose(arg))/2
      elif isinstance(arg,escript.Data):
-       return escript_nonsymmetric(arg)
+       if arg.getRank()==2:
+         if not (arg.getShape()[0]==arg.getShape()[1]):
+            raise ValueError,"argument must be square."
+         return arg._nonsymmetric()
+       elif arg.getRank()==4:
+         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
+            raise ValueError,"argument must be square."
+         return arg._nonsymmetric()
+       else:
+         raise ValueError,"rank 2 or 4 is required."
      elif isinstance(arg,float):
        return arg
      elif isinstance(arg,int):
-Line 3374 
 def nonsymmetric(arg):
+Line 3364 
 def nonsymmetric(arg):
      else:
        raise TypeError,"nonsymmetric: Unknown argument type."
- 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
  def inverse(arg):
      """
      returns the inverse of the square matrix arg.
      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
-     @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
+     @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
-     @remark: for L{escript.Data} objects the dimension is restricted to 3.
+     @note: for L{escript.Data} objects the dimension is restricted to 3.
      """
+     import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
      if isinstance(arg,numarray.NumArray):
        return numarray.linear_algebra.inverse(arg)
      elif isinstance(arg,escript.Data):
-Line 3498 
 class Inverse_Symbol(DependendSymbol):
+Line 3467 
 class Inverse_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str"  or format=="text":
           return "inverse(%s)"%argstrs[0]
-Line 3538 
 class Inverse_Symbol(DependendSymbol):
+Line 3507 
 class Inverse_Symbol(DependendSymbol):
        if arg==self:
           return identity(self.getShape())
        else:
-          return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)
+          return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
  def eigenvalues(arg):
      """
-Line 3549 
 def eigenvalues(arg):
+Line 3518 
 def eigenvalues(arg):
      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
      @return: the eigenvalues in increasing order.
      @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
-     @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.
      """
      if isinstance(arg,numarray.NumArray):
        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
-Line 3610 
 def eigenvalues(arg):
+Line 3579 
 def eigenvalues(arg):
  def eigenvalues_and_eigenvectors(arg):
      """
-     returns the eigenvalues of the square matrix arg.
+     returns the eigenvalues and eigenvectors of the square matrix arg.
      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
                  arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
-     @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @type arg: L{escript.Data}
-     @return: the eigenvalues in increasing order.
+     @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
-     @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
-     @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
+              the eigenvector coresponding to the i-th eigenvalue.
+     @rtype: L{tuple} of L{escript.Data}.
+     @note: The dimension is restricted to 3.
      """
      if isinstance(arg,numarray.NumArray):
        raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
-Line 3690 
 class Add_Symbol(DependendSymbol):
+Line 3661 
 class Add_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="str" or format=="text":
           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
-Line 3789 
 class Mult_Symbol(DependendSymbol):
+Line 3760 
 class Mult_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="str" or format=="text":
           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
-Line 3894 
 class Quotient_Symbol(DependendSymbol):
+Line 3865 
 class Quotient_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="str" or format=="text":
           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
-Line 3998 
 class Power_Symbol(DependendSymbol):
+Line 3969 
 class Power_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str" or format=="text":
           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
-Line 4087 
 def clip(arg,minval=0.,maxval=1.):
+Line 4058 
 def clip(arg,minval=0.,maxval=1.):
      @param arg: argument
      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
      @param minval: lower range
-     @type arg: C{float}
+     @type minval: C{float}
      @param maxval: upper range
-     @type arg: C{float}
+     @type maxval: C{float}
      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
               less then maxval are unchanged.
      @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 4085 
 def inner(arg0,arg1):
      @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 inner product of arg0 and arg1 at each data point
+     @return: the inner product of arg0 and arg1 at each data point
      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input
      @raise ValueError: if the shapes of the arguments are not identical
      """
-Line 4124 
 def inner(arg0,arg1):
+Line 4095 
 def inner(arg0,arg1):
          raise ValueError,"inner: shape of arguments does not match"
      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
+ def outer(arg0,arg1):
+     """
+     the outer product of the two argument:
+     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
+     """
+     return generalTensorProduct(arg0,arg1,axis_offset=0)
  def matrixmult(arg0,arg1):
      """
+     see L{matrix_mult}
+     """
+     return matrix_mult(arg0,arg1)
+ def matrix_mult(arg0,arg1):
+     """
      matrix-matrix or matrix-vector product of the two argument:
      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
-           or
+     or
      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
-     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.
      @param arg0: first argument of rank 2
      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
-Line 4152 
 def matrixmult(arg0,arg1):
+Line 4149 
 def matrixmult(arg0,arg1):
          raise ValueError,"second argument must have rank 1 or 2"
      return generalTensorProduct(arg0,arg1,axis_offset=1)
- def outer(arg0,arg1):
+ def tensormult(arg0,arg1):
      """
-     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
      """
-     return generalTensorProduct(arg0,arg1,axis_offset=0)
+     return tensor_mult(arg0,arg1)
+ def tensor_mult(arg0,arg1):
- def tensormult(arg0,arg1):
      """
      the tensor product of the two argument:
      for arg0 of rank 2 this is
      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
-                  or
+     or
      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
-Line 4189 
 def tensormult(arg0,arg1):
+Line 4172 
 def tensormult(arg0,arg1):
      out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2,s3]
-                  or
+     or
      out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2]
-                  or
+     or
      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]
-     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
-     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.
      @param arg0: first argument of rank 2 or 4
      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
-Line 4214 
 def tensormult(arg0,arg1):
+Line 4197 
 def tensormult(arg0,arg1):
      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
         return generalTensorProduct(arg0,arg1,axis_offset=2)
      else:
-         raise ValueError,"tensormult: first argument must have rank 2 or 4"
+         raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
  def generalTensorProduct(arg0,arg1,axis_offset=0):
      """
-Line 4222 
 def generalTensorProduct(arg0,arg1,axis_
+Line 4205 
 def generalTensorProduct(arg0,arg1,axis_
      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
-     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:]
-     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]
-     in the second case the two last dimensions of arg0 must match the shape of arg1.
+         - r runs trough arg0.Shape[:axis_offset]
+         - t runs through arg1.Shape[axis_offset:]
      @param arg0: first argument
      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
-     @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
+     @param arg1: second argument
      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
      @return: the general tensor product of arg0 and arg1 at each data point.
      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
-Line 4244 
 def generalTensorProduct(arg0,arg1,axis_
+Line 4226 
 def generalTensorProduct(arg0,arg1,axis_
             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
         else:
             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
-                raise ValueError,"generalTensorProduct: dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
+                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)
             arg0_c=arg0.copy()
             arg1_c=arg1.copy()
             sh0,sh1=arg0.shape,arg1.shape
-Line 4270 
 def generalTensorProduct(arg0,arg1,axis_
+Line 4252 
 def generalTensorProduct(arg0,arg1,axis_
  class GeneralTensorProduct_Symbol(DependendSymbol):
     """
-    Symbol representing the quotient of two arguments.
+    Symbol representing the general tensor product of two arguments
     """
     def __init__(self,arg0,arg1,axis_offset=0):
         """
-        initialization of L{Symbol} representing the quotient of two arguments
+        initialization of L{Symbol} representing the general tensor product of two arguments.
-        @param arg0: numerator
+        @param arg0: first argument
         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @param arg1: denominator
+        @param arg1: second argument
         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @raise ValueError: if both arguments do not have the same shape.
+        @raise ValueError: illegal dimension
         @note: if both arguments have a spatial dimension, they must equal.
         """
         sh_arg0=pokeShape(arg0)
-Line 4303 
 class GeneralTensorProduct_Symbol(Depend
+Line 4285 
 class GeneralTensorProduct_Symbol(Depend
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str" or format=="text":
           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
-Line 4331 
 class GeneralTensorProduct_Symbol(Depend
+Line 4313 
 class GeneralTensorProduct_Symbol(Depend
           args=self.getSubstitutedArguments(argvals)
           return generalTensorProduct(args[0],args[1],args[2])
+ def escript_generalTensorProductNew(arg0,arg1,axis_offset):
+     "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
+     # calculate the return shape:
+     shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
+     shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
+     shape10=arg1.getShape()[:axis_offset]
+     shape1=arg1.getShape()[axis_offset:]
+     if not shape01==shape10:
+         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)
+     # Figure out which functionspace to use (look at where operator+ is defined maybe in BinaryOp.h to get the logic for this)
+     # fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
+     out=GeneralTensorProduct(arg0, arg1, axis_offset)
+     return out
  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct
      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
      # calculate the return shape:
-Line 4374 
 def escript_generalTensorProduct(arg0,ar
+Line 4370 
 def escript_generalTensorProduct(arg0,ar
      return out
+ def transposed_matrix_mult(arg0,arg1):
+     """
+     transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
+     out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
+     or
+     out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
+     The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
+     The first dimension of arg0 and arg1 must match.
+     @param arg0: first argument of rank 2
+     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @param arg1: second argument of at least rank 1
+     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @return: the product of the transposed of arg0 and arg1 at each data point
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     @raise ValueError: if the shapes of the arguments are not appropriate
+     """
+     sh0=pokeShape(arg0)
+     sh1=pokeShape(arg1)
+     if not len(sh0)==2 :
+         raise ValueError,"first argument must have rank 2"
+     if not len(sh1)==2 and not len(sh1)==1:
+         raise ValueError,"second argument must have rank 1 or 2"
+     return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
+ def transposed_tensor_mult(arg0,arg1):
+     """
+     the tensor product of the transposed of the first and the second argument
+     for arg0 of rank 2 this is
+     out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
+     or
+     out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
+     and for arg0 of rank 4 this is
+     out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
+     or
+     out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
+     or
+     out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
+     In the first case the the first dimension of arg0 and the first dimension of arg1 must match and
+     in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
+     The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
+     @param arg0: first argument of rank 2 or 4
+     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
+     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     """
+     sh0=pokeShape(arg0)
+     sh1=pokeShape(arg1)
+     if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
+        return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
+     elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
+        return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
+     else:
+         raise ValueError,"first argument must have rank 2 or 4"
+ def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
+     """
+     generalized tensor product of transposed of arg0 and arg1:
+     out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
+     where
+         - s runs through arg0.Shape[axis_offset:]
+         - r runs trough arg0.Shape[:axis_offset]
+         - t runs through arg1.Shape[axis_offset:]
+     The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
+     to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
+     @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 general tensor product of transposed(arg0) and arg1 at each data point.
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     """
+     if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
+     arg0,arg1=matchType(arg0,arg1)
+     # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
+     if isinstance(arg0,numarray.NumArray):
+        if isinstance(arg1,Symbol):
+            return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
+        else:
+            if not arg0.shape[:axis_offset]==arg1.shape[: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)
+            arg0_c=arg0.copy()
+            arg1_c=arg1.copy()
+            sh0,sh1=arg0.shape,arg1.shape
+            d0,d1,d01=1,1,1
+            for i in sh0[axis_offset:]: d0*=i
+            for i in sh1[axis_offset:]: d1*=i
+            for i in sh0[:axis_offset]: d01*=i
+            arg0_c.resize((d01,d0))
+            arg1_c.resize((d01,d1))
+            out=numarray.zeros((d0,d1),numarray.Float64)
+            for i0 in range(d0):
+                     for i1 in range(d1):
+                          out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
+            out.resize(sh0[axis_offset:]+sh1[axis_offset:])
+            return out
+     elif isinstance(arg0,escript.Data):
+        if isinstance(arg1,Symbol):
+            return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
+        else:
+            return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
+     else:
+        return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
+ class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
+    """
+    Symbol representing the general tensor product of the transposed of arg0 and arg1
+    """
+    def __init__(self,arg0,arg1,axis_offset=0):
+        """
+        initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
+        @param arg0: first argument
+        @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
+        @param arg1: second argument
+        @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
+        @raise ValueError: inconsistent dimensions of arguments.
+        @note: if both arguments have a spatial dimension, they must equal.
+        """
+        sh_arg0=pokeShape(arg0)
+        sh_arg1=pokeShape(arg1)
+        sh01=sh_arg0[:axis_offset]
+        sh10=sh_arg1[:axis_offset]
+        sh0=sh_arg0[axis_offset:]
+        sh1=sh_arg1[axis_offset:]
+        if not sh01==sh10:
+            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)
+        DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
+    def getMyCode(self,argstrs,format="escript"):
+       """
+       returns a program code that can be used to evaluate the symbol.
+       @param argstrs: gives for each argument a string representing the argument for the evaluation.
+       @type argstrs: C{list} of length 2 of C{str}.
+       @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
+       @type format: C{str}
+       @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
+       @rtype: C{str}
+       @raise NotImplementedError: if the requested format is not available
+       """
+       if format=="escript" or format=="str" or format=="text":
+          return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
+       else:
+          raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
+    def substitute(self,argvals):
+       """
+       assigns new values to symbols in the definition of the symbol.
+       The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
+       @param argvals: new values assigned to symbols
+       @type argvals: C{dict} with keywords of type L{Symbol}.
+       @return: result of the substitution process. Operations are executed as much as possible.
+       @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
+       @raise TypeError: if a value for a L{Symbol} cannot be substituted.
+       """
+       if argvals.has_key(self):
+          arg=argvals[self]
+          if self.isAppropriateValue(arg):
+             return arg
+          else:
+             raise TypeError,"%s: new value is not appropriate."%str(self)
+       else:
+          args=self.getSubstitutedArguments(argvals)
+          return generalTransposedTensorProduct(args[0],args[1],args[2])
+ def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
+     "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
+     # calculate the return shape:
+     shape01=arg0.getShape()[:axis_offset]
+     shape10=arg1.getShape()[:axis_offset]
+     shape0=arg0.getShape()[axis_offset:]
+     shape1=arg1.getShape()[axis_offset:]
+     if not shape01==shape10:
+         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)
+     # whatr function space should be used? (this here is not good!)
+     fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
+     # create return value:
+     out=escript.Data(0.,tuple(shape0+shape1),fs)
+     #
+     s0=[[]]
+     for k in shape0:
+           s=[]
+           for j in s0:
+                 for i in range(k): s.append(j+[slice(i,i)])
+           s0=s
+     s1=[[]]
+     for k in shape1:
+           s=[]
+           for j in s1:
+                 for i in range(k): s.append(j+[slice(i,i)])
+           s1=s
+     s01=[[]]
+     for k in shape01:
+           s=[]
+           for j in s01:
+                 for i in range(k): s.append(j+[slice(i,i)])
+           s01=s
+     for i0 in s0:
+        for i1 in s1:
+          s=escript.Scalar(0.,fs)
+          for i01 in s01:
+             s+=arg0.__getitem__(tuple(i01+i0))*arg1.__getitem__(tuple(i01+i1))
+          out.__setitem__(tuple(i0+i1),s)
+     return out
+ def matrix_transposed_mult(arg0,arg1):
+     """
+     matrix-transposed(matrix) product of the two argument:
+     out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
+     The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
+     The last dimensions of arg0 and arg1 must match.
+     @param arg0: first argument of rank 2
+     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @param arg1: second argument of rank 2
+     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @return: the product of arg0 and the transposed of arg1 at each data point
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     @raise ValueError: if the shapes of the arguments are not appropriate
+     """
+     sh0=pokeShape(arg0)
+     sh1=pokeShape(arg1)
+     if not len(sh0)==2 :
+         raise ValueError,"first argument must have rank 2"
+     if not len(sh1)==2 and not len(sh1)==1:
+         raise ValueError,"second argument must have rank 1 or 2"
+     return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
+ def tensor_transposed_mult(arg0,arg1):
+     """
+     the tensor product of the first and the transpose of the second argument
+     for arg0 of rank 2 this is
+     out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
+     and for arg0 of rank 4 this is
+     out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
+     or
+     out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
+     In the first case the the second dimension of arg0 and arg1 must match and
+     in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
+     The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
+     @param arg0: first argument of rank 2 or 4
+     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
+     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     """
+     sh0=pokeShape(arg0)
+     sh1=pokeShape(arg1)
+     if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
+        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
+     elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
+        return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
+     else:
+         raise ValueError,"first argument must have rank 2 or 4"
+ def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
+     """
+     generalized tensor product of transposed of arg0 and arg1:
+     out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
+     where
+         - s runs through arg0.Shape[:arg0.Rank-axis_offset]
+         - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
+         - t runs through arg1.Shape[arg1.Rank-axis_offset:]
+     The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
+     to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
+     @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 general tensor product of transposed(arg0) and arg1 at each data point.
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     """
+     if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
+     arg0,arg1=matchType(arg0,arg1)
+     # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
+     if isinstance(arg0,numarray.NumArray):
+        if isinstance(arg1,Symbol):
+            return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
+        else:
+            if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-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)
+            arg0_c=arg0.copy()
+            arg1_c=arg1.copy()
+            sh0,sh1=arg0.shape,arg1.shape
+            d0,d1,d01=1,1,1
+            for i in sh0[:arg0.rank-axis_offset]: d0*=i
+            for i in sh1[:arg1.rank-axis_offset]: d1*=i
+            for i in sh1[arg1.rank-axis_offset:]: d01*=i
+            arg0_c.resize((d0,d01))
+            arg1_c.resize((d1,d01))
+            out=numarray.zeros((d0,d1),numarray.Float64)
+            for i0 in range(d0):
+                     for i1 in range(d1):
+                          out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
+            out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
+            return out
+     elif isinstance(arg0,escript.Data):
+        if isinstance(arg1,Symbol):
+            return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
+        else:
+            return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
+     else:
+        return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
+ class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
+    """
+    Symbol representing the general tensor product of arg0 and the transpose of arg1
+    """
+    def __init__(self,arg0,arg1,axis_offset=0):
+        """
+        initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
+        @param arg0: first argument
+        @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
+        @param arg1: second argument
+        @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
+        @raise ValueError: inconsistent dimensions of arguments.
+        @note: if both arguments have a spatial dimension, they must equal.
+        """
+        sh_arg0=pokeShape(arg0)
+        sh_arg1=pokeShape(arg1)
+        sh0=sh_arg0[:len(sh_arg0)-axis_offset]
+        sh01=sh_arg0[len(sh_arg0)-axis_offset:]
+        sh10=sh_arg1[len(sh_arg1)-axis_offset:]
+        sh1=sh_arg1[:len(sh_arg1)-axis_offset]
+        if not sh01==sh10:
+            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)
+        DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
+    def getMyCode(self,argstrs,format="escript"):
+       """
+       returns a program code that can be used to evaluate the symbol.
+       @param argstrs: gives for each argument a string representing the argument for the evaluation.
+       @type argstrs: C{list} of length 2 of C{str}.
+       @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
+       @type format: C{str}
+       @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
+       @rtype: C{str}
+       @raise NotImplementedError: if the requested format is not available
+       """
+       if format=="escript" or format=="str" or format=="text":
+          return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
+       else:
+          raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
+    def substitute(self,argvals):
+       """
+       assigns new values to symbols in the definition of the symbol.
+       The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
+       @param argvals: new values assigned to symbols
+       @type argvals: C{dict} with keywords of type L{Symbol}.
+       @return: result of the substitution process. Operations are executed as much as possible.
+       @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
+       @raise TypeError: if a value for a L{Symbol} cannot be substituted.
+       """
+       if argvals.has_key(self):
+          arg=argvals[self]
+          if self.isAppropriateValue(arg):
+             return arg
+          else:
+             raise TypeError,"%s: new value is not appropriate."%str(self)
+       else:
+          args=self.getSubstitutedArguments(argvals)
+          return generalTensorTransposedProduct(args[0],args[1],args[2])
+ def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
+     "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
+     # calculate the return shape:
+     shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
+     shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
+     shape10=arg1.getShape()[arg1.getRank()-axis_offset:]
+     shape1=arg1.getShape()[:arg1.getRank()-axis_offset]
+     if not shape01==shape10:
+         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)
+     # whatr function space should be used? (this here is not good!)
+     fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
+     # create return value:
+     out=escript.Data(0.,tuple(shape0+shape1),fs)
+     #
+     s0=[[]]
+     for k in shape0:
+           s=[]
+           for j in s0:
+                 for i in range(k): s.append(j+[slice(i,i)])
+           s0=s
+     s1=[[]]
+     for k in shape1:
+           s=[]
+           for j in s1:
+                 for i in range(k): s.append(j+[slice(i,i)])
+           s1=s
+     s01=[[]]
+     for k in shape01:
+           s=[]
+           for j in s01:
+                 for i in range(k): s.append(j+[slice(i,i)])
+           s01=s
+     for i0 in s0:
+        for i1 in s1:
+          s=escript.Scalar(0.,fs)
+          for i01 in s01:
+             s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i1+i01))
+          out.__setitem__(tuple(i0+i1),s)
+     return out
  #=========================================================
  #  functions dealing with spatial dependency
  #=========================================================
-Line 4394 
 def grad(arg,where=None):
+Line 4849 
 def grad(arg,where=None):
                    If not present or C{None} an appropriate default is used.
      @type where: C{None} or L{escript.FunctionSpace}
      @return: gradient of arg.
-     @rtype:  L{escript.Data} or L{Symbol}
+     @rtype: L{escript.Data} or L{Symbol}
      """
      if isinstance(arg,Symbol):
         return Grad_Symbol(arg,where)
-Line 4434 
 class Grad_Symbol(DependendSymbol):
+Line 4889 
 class Grad_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str"  or format=="text":
           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
-Line 4487 
 def integrate(arg,where=None):
+Line 4942 
 def integrate(arg,where=None):
                    If not present or C{None} an appropriate default is used.
      @type where: C{None} or L{escript.FunctionSpace}
      @return: integral of arg.
-     @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}
+     @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
      """
      if isinstance(arg,Symbol):
         return Integrate_Symbol(arg,where)
-Line 4525 
 class Integrate_Symbol(DependendSymbol):
+Line 4980 
 class Integrate_Symbol(DependendSymbol):
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str"  or format=="text":
           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
-Line 4577 
 def interpolate(arg,where):
+Line 5032 
 def interpolate(arg,where):
      @param where: FunctionSpace to be interpolated to
      @type where: L{escript.FunctionSpace}
      @return: interpolated argument
-     @rtype:  C{escript.Data} or L{Symbol}
+     @rtype: C{escript.Data} or L{Symbol}
      """
      if isinstance(arg,Symbol):
         return Interpolate_Symbol(arg,where)
-Line 4608 
 class Interpolate_Symbol(DependendSymbol
+Line 5063 
 class Interpolate_Symbol(DependendSymbol
        @type format: C{str}
        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
        @rtype: C{str}
-       @raise: NotImplementedError: if the requested format is not available
+       @raise NotImplementedError: if the requested format is not available
        """
        if format=="escript" or format=="str"  or format=="text":
           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
-Line 4661 
 def div(arg,where=None):
+Line 5116 
 def div(arg,where=None):
                    If not present or C{None} an appropriate default is used.
      @type where: C{None} or L{escript.FunctionSpace}
      @return: divergence of arg.
-     @rtype:  L{escript.Data} or L{Symbol}
+     @rtype: L{escript.Data} or L{Symbol}
      """
      if isinstance(arg,Symbol):
          dim=arg.getDim()
-Line 4683 
 def jump(arg,domain=None):
+Line 5138 
 def jump(arg,domain=None):
                     the domain of arg is used. If arg is a L{Symbol} the domain must be present.
      @type domain: C{None} or L{escript.Domain}
      @return: jump of arg
-     @rtype:  L{escript.Data} or L{Symbol}
+     @rtype: L{escript.Data} or L{Symbol}
      """
      if domain==None: domain=arg.getDomain()
      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
-Line 4695 
 def L2(arg):
+Line 5150 
 def L2(arg):
      @param arg: function which L2 to be calculated.
      @type arg: L{escript.Data} or L{Symbol}
      @return: L2 norm of arg.
-     @rtype:  L{float} or L{Symbol}
+     @rtype: L{float} or L{Symbol}
      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
      """
      return sqrt(integrate(inner(arg,arg)))

 Legend:



Removed from v.587
 


changed lines


 
Added in v.804
 Legend:



Removed from v.587
 


changed lines


 
Added in v.804
-Removed from v.587
+Added in v.804

	ViewVC Help
Powered by ViewVC 1.1.26