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

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

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

revision 341 by gross, Mon Dec 12 05:26:10 2005 UTC revision 804 by gross, Thu Aug 10 01:12:16 2006 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2    
3  """  """
4  Utility functions for escript  Utility functions for escript
5    
 @remark:  This module is under construction and is still tested!!!  
   
6  @var __author__: name of author  @var __author__: name of author
7  @var __licence__: licence agreement  @var __copyright__: copyrights
8    @var __license__: licence agreement
9  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
10  @var __version__: version  @var __version__: version
11  @var __date__: date of the version  @var __date__: date of the version
12  """  """
13                                                                                                                                                                                                                                                                                                                                                                                                            
14  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
15  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                        http://www.access.edu.au
17                    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19                 http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
21  __version__="$Revision: 329 $"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
23    
24    
# Line 33  import numarray Line 27  import numarray
27  import escript  import escript
28  import os  import os
29    
 # missing tests:  
   
 # def pokeShape(arg):  
 # def pokeDim(arg):  
 # def commonShape(arg0,arg1):  
 # def commonDim(*args):  
 # def testForZero(arg):  
 # def matchType(arg0=0.,arg1=0.):  
 # def matchShape(arg0,arg1):  
   
 # def maximum(arg0,arg1):  
 # def minimum(arg0,arg1):  
   
 # def transpose(arg,axis=None):  
 # def trace(arg,axis0=0,axis1=1):  
 # def reorderComponents(arg,index):  
   
 # def integrate(arg,where=None):  
 # def interpolate(arg,where):  
 # def div(arg,where=None):  
 # def grad(arg,where=None):  
   
 #  
 # slicing: get  
 #          set  
 #  
 # and derivatives  
   
30  #=========================================================  #=========================================================
31  #   some helpers:  #   some helpers:
32  #=========================================================  #=========================================================
# Line 68  def saveVTK(filename,domain=None,**data) Line 34  def saveVTK(filename,domain=None,**data)
34      """      """
35      writes a L{Data} objects into a files using the the VTK XML file format.      writes a L{Data} objects into a files using the the VTK XML file format.
36    
37      Example:      Example::
38    
39         tmp=Scalar(..)         tmp=Scalar(..)
40         v=Vector(..)         v=Vector(..)
# Line 96  def saveDX(filename,domain=None,**data): Line 62  def saveDX(filename,domain=None,**data):
62      """      """
63      writes a L{Data} objects into a files using the the DX file format.      writes a L{Data} objects into a files using the the DX file format.
64    
65      Example:      Example::
66    
67         tmp=Scalar(..)         tmp=Scalar(..)
68         v=Vector(..)         v=Vector(..)
# Line 125  def kronecker(d=3): Line 91  def kronecker(d=3):
91     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
92    
93     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
94     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
95     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
96     @rtype d: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2.
    @remark: the function is identical L{identity}  
97     """     """
98     return identityTensor(d)     return identityTensor(d)
99    
# Line 143  def identity(shape=()): Line 108  def identity(shape=()):
108     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
109     """     """
110     if len(shape)>0:     if len(shape)>0:
111        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
112        if len(shape)==1:        if len(shape)==1:
113            for i0 in range(shape[0]):            for i0 in range(shape[0]):
114               out[i0,i0]=1.               out[i0,i0]=1.
   
115        elif len(shape)==2:        elif len(shape)==2:
116            for i0 in range(shape[0]):            for i0 in range(shape[0]):
117               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 163  def identityTensor(d=3): Line 127  def identityTensor(d=3):
127     return the dxd identity matrix     return the dxd identity matrix
128    
129     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
130     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
131     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
132     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
133     """     """
134     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
135        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
136     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
137           return identity((d.getDim(),))
138       else:
139           return identity((d,))
140    
141  def identityTensor4(d=3):  def identityTensor4(d=3):
142     """     """
# Line 178  def identityTensor4(d=3): Line 145  def identityTensor4(d=3):
145     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
146     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
147     @return: the object u of rank 4 with M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise     @return: the object u of rank 4 with M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise
148     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
149     """     """
150     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
151        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
152     return identity((d,d))     elif isinstance(d,escript.Domain):
153           return identity((d.getDim(),d.getDim()))
154       else:
155           return identity((d,d))
156    
157  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
158     """     """
# Line 191  def unitVector(i=0,d=3): Line 161  def unitVector(i=0,d=3):
161     @param i: index     @param i: index
162     @type i: C{int}     @type i: C{int}
163     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
164     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
165     @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise     @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise
166     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
167     """     """
168     return kronecker(d)[i]     return kronecker(d)[i]
169    
# Line 249  def inf(arg): Line 219  def inf(arg):
219    
220      @param arg: argument      @param arg: argument
221      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
222      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
223      @rtype: C{float}      @rtype: C{float}
224      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
225      """      """
# Line 338  def commonDim(*args): Line 308  def commonDim(*args):
308      """      """
309      identifies, if possible, the spatial dimension across a set of objects which may or my not have a spatial dimension.      identifies, if possible, the spatial dimension across a set of objects which may or my not have a spatial dimension.
310    
311      @param *args: given objects      @param args: given objects
312      @return: the spatial dimension of the objects with identifiable dimension (see L{pokeDim}). If none the objects has      @return: the spatial dimension of the objects with identifiable dimension (see L{pokeDim}). If none the objects has
313               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
314      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 360  def testForZero(arg): Line 330  def testForZero(arg):
330    
331      @param arg: a given object      @param arg: a given object
332      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}
333      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
334      @rtype : C{bool}      @rtype: C{bool}
335      """      """
336      try:      if isinstance(arg,numarray.NumArray):
337           return not Lsup(arg)>0.
338        elif isinstance(arg,escript.Data):
339           return False
340        elif isinstance(arg,float):
341           return not Lsup(arg)>0.
342        elif isinstance(arg,int):
343         return not Lsup(arg)>0.         return not Lsup(arg)>0.
344      except TypeError:      elif isinstance(arg,Symbol):
345           return False
346        else:
347         return False         return False
348    
349  def matchType(arg0=0.,arg1=0.):  def matchType(arg0=0.,arg1=0.):
# Line 386  def matchType(arg0=0.,arg1=0.): Line 364  def matchType(arg0=0.,arg1=0.):
364         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
365            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
366         elif isinstance(arg1,float):         elif isinstance(arg1,float):
367            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
368         elif isinstance(arg1,int):         elif isinstance(arg1,int):
369            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
370         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
371            pass            pass
372         else:         else:
# Line 412  def matchType(arg0=0.,arg1=0.): Line 390  def matchType(arg0=0.,arg1=0.):
390         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
391            pass            pass
392         elif isinstance(arg1,float):         elif isinstance(arg1,float):
393            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
394         elif isinstance(arg1,int):         elif isinstance(arg1,int):
395            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
396         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
397            pass            pass
398         else:         else:
399            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
400      elif isinstance(arg0,float):      elif isinstance(arg0,float):
401         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
402            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
403         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
404            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
405         elif isinstance(arg1,float):         elif isinstance(arg1,float):
406            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
407            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
408         elif isinstance(arg1,int):         elif isinstance(arg1,int):
409            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
410            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
411         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
412            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
413         else:         else:
414            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
415      elif isinstance(arg0,int):      elif isinstance(arg0,int):
416         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
417            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
418         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
419            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
420         elif isinstance(arg1,float):         elif isinstance(arg1,float):
421            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
422            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
423         elif isinstance(arg1,int):         elif isinstance(arg1,int):
424            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
425            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
426         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
427            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
428         else:         else:
429            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
430      else:      else:
# Line 456  def matchType(arg0=0.,arg1=0.): Line 434  def matchType(arg0=0.,arg1=0.):
434    
435  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
436      """      """
437            return representations of arg0 amd arg1 which ahve the same shape
   
     If shape is not given the shape "largest" shape of args is used.  
438    
439      @param args: a given ob      @param arg0: a given object
440      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
441      @return: True if the argument is identical to zero.      @param arg1: a given object
442      @rtype: C{list} of C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
443        @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
444        @rtype: C{tuple}
445      """      """
446      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
447      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
448      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
449      if len(sh0)<len(sh):      if len(sh0)<len(sh):
450         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
451      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
452         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
453      else:      else:
454         return arg0,arg1         return arg0,arg1
455  #=========================================================  #=========================================================
# Line 491  class Symbol(object): Line 469  class Symbol(object):
469         Creates an instance of a symbol of a given shape. The symbol may depending on a list of arguments args which may be         Creates an instance of a symbol of a given shape. The symbol may depending on a list of arguments args which may be
470         symbols or any other object.         symbols or any other object.
471    
472         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
473         @type arg: C{list}         @type args: C{list}
474         @param shape: the shape         @param shape: the shape
475         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
476         @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
# Line 535  class Symbol(object): Line 513  class Symbol(object):
513         """         """
514         the shape of the symbol.         the shape of the symbol.
515    
516         @return : the shape of the symbol.         @return: the shape of the symbol.
517         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
518         """         """
519         return self.__shape         return self.__shape
# Line 544  class Symbol(object): Line 522  class Symbol(object):
522         """         """
523         the spatial dimension         the spatial dimension
524    
525         @return : the spatial dimension         @return: the spatial dimension
526         @rtype: C{int} if the dimension is defined. Otherwise C{None} is returned.         @rtype: C{int} if the dimension is defined. Otherwise C{None} is returned.
527         """         """
528         return self.__dim         return self.__dim
# Line 568  class Symbol(object): Line 546  class Symbol(object):
546         """         """
547         substitutes symbols in the arguments of this object and returns the result as a list.         substitutes symbols in the arguments of this object and returns the result as a list.
548    
549         @param argvals: L{Symbols} and their substitutes. The L{Symbol} u in the expression defining this object is replaced by argvals[u].         @param argvals: L{Symbol} and their substitutes. The L{Symbol} u in the expression defining this object is replaced by argvals[u].
550         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
551         @rtype: C{list} of objects         @rtype: C{list} of objects
552         @return: list of the object assigned to the arguments through substitution or for the arguments which are not L{Symbols} the value assigned to the argument at instantiation.         @return: list of the object assigned to the arguments through substitution or for the arguments which are not L{Symbol} the value assigned to the argument at instantiation.
553         """         """
554         out=[]         out=[]
555         for a in self.getArgument():         for a in self.getArgument():
# Line 597  class Symbol(object): Line 575  class Symbol(object):
575            else:            else:
576                s=pokeShape(s)+arg.getShape()                s=pokeShape(s)+arg.getShape()
577                if len(s)>0:                if len(s)>0:
578                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
579                else:                else:
580                   out.append(a)                   out.append(a)
581         return out         return out
# Line 687  class Symbol(object): Line 665  class Symbol(object):
665         else:         else:
666            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
667            if len(s)>0:            if len(s)>0:
668               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
669            else:            else:
670               return 0.               return 0.
671    
# Line 695  class Symbol(object): Line 673  class Symbol(object):
673         """         """
674         returns -self.         returns -self.
675    
676         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
677         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
678         """         """
679         return self*(-1.)         return self*(-1.)
# Line 704  class Symbol(object): Line 682  class Symbol(object):
682         """         """
683         returns +self.         returns +self.
684    
685         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
686         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
687         """         """
688         return self*(1.)         return self*(1.)
689    
690     def __abs__(self):     def __abs__(self):
691         """         """
692         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
693         """         """
694         return Abs_Symbol(self)         return Abs_Symbol(self)
695    
# Line 721  class Symbol(object): Line 699  class Symbol(object):
699    
700         @param other: object to be added to this object         @param other: object to be added to this object
701         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
702         @return:  a S{Symbol} representing the sum of this object and C{other}         @return:  a L{Symbol} representing the sum of this object and C{other}
703         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
704         """         """
705         return add(self,other)         return add(self,other)
# Line 732  class Symbol(object): Line 710  class Symbol(object):
710    
711         @param other: object this object is added to         @param other: object this object is added to
712         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
713         @return: a S{Symbol} representing the sum of C{other} and this object object         @return: a L{Symbol} representing the sum of C{other} and this object object
714         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
715         """         """
716         return add(other,self)         return add(other,self)
# Line 743  class Symbol(object): Line 721  class Symbol(object):
721    
722         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
723         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
724         @return: a S{Symbol} representing the difference of C{other} and this object         @return: a L{Symbol} representing the difference of C{other} and this object
725         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
726         """         """
727         return add(self,-other)         return add(self,-other)
# Line 754  class Symbol(object): Line 732  class Symbol(object):
732    
733         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
734         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
735         @return: a S{Symbol} representing the difference of this object and C{other}.         @return: a L{Symbol} representing the difference of this object and C{other}.
736         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
737         """         """
738         return add(-self,other)         return add(-self,other)
# Line 765  class Symbol(object): Line 743  class Symbol(object):
743    
744         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
745         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
746         @return: a S{Symbol} representing the product of the object and C{other}.         @return: a L{Symbol} representing the product of the object and C{other}.
747         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
748         """         """
749         return mult(self,other)         return mult(self,other)
# Line 776  class Symbol(object): Line 754  class Symbol(object):
754    
755         @param other: object this object is multiplied with         @param other: object this object is multiplied with
756         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
757         @return: a S{Symbol} representing the product of C{other} and the object.         @return: a L{Symbol} representing the product of C{other} and the object.
758         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
759         """         """
760         return mult(other,self)         return mult(other,self)
# Line 787  class Symbol(object): Line 765  class Symbol(object):
765    
766         @param other: object dividing this object         @param other: object dividing this object
767         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
768         @return: a S{Symbol} representing the quotient of this object and C{other}         @return: a L{Symbol} representing the quotient of this object and C{other}
769         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
770         """         """
771         return quotient(self,other)         return quotient(self,other)
# Line 798  class Symbol(object): Line 776  class Symbol(object):
776    
777         @param other: object dividing this object         @param other: object dividing this object
778         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
779         @return: a S{Symbol} representing the quotient of C{other} and this object         @return: a L{Symbol} representing the quotient of C{other} and this object
780         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
781         """         """
782         return quotient(other,self)         return quotient(other,self)
# Line 809  class Symbol(object): Line 787  class Symbol(object):
787    
788         @param other: exponent         @param other: exponent
789         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
790         @return: a S{Symbol} representing the power of this object to C{other}         @return: a L{Symbol} representing the power of this object to C{other}
791         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
792         """         """
793         return power(self,other)         return power(self,other)
# Line 820  class Symbol(object): Line 798  class Symbol(object):
798    
799         @param other: basis         @param other: basis
800         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
801         @return: a S{Symbol} representing the power of C{other} to this object         @return: a L{Symbol} representing the power of C{other} to this object
802         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
803         """         """
804         return power(other,self)         return power(other,self)
805    
806       def __getitem__(self,index):
807           """
808           returns the slice defined by index
809    
810           @param index: defines a
811           @type index: C{slice} or C{int} or a C{tuple} of them
812           @return: a L{Symbol} representing the slice defined by index
813           @rtype: L{DependendSymbol}
814           """
815           return GetSlice_Symbol(self,index)
816    
817  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
818     """     """
819     DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.     DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.
820     Two DependendSymbol are equal if they have the same shape, the same arguments and one of them has an unspecified spatial dimension or the spatial dimension is identical       Two DependendSymbol are equal if they have the same shape, the same arguments and one of them has an unspecified spatial dimension or the spatial dimension is identical  
821        
822     Example:     Example::
823        
824     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
825     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
826     print u1==u2       print u1==u2
827     False       False
828        
829        but       but::
830    
831     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
832     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
833     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
834     print u1==u2, u1==u3       print u1==u2, u1==u3
835     True False       True False
836    
837     @note: DependendSymbol should be used as return value of functions with L{Symbol} arguments. This will allow the optimizer to remove redundant function calls.     @note: DependendSymbol should be used as return value of functions with L{Symbol} arguments. This will allow the optimizer to remove redundant function calls.
838     """     """
# Line 875  class DependendSymbol(Symbol): Line 864  class DependendSymbol(Symbol):
864  #=========================================================  #=========================================================
865  #  Unary operations prserving the shape  #  Unary operations prserving the shape
866  #========================================================  #========================================================
867    class GetSlice_Symbol(DependendSymbol):
868       """
869       L{Symbol} representing getting a slice for a L{Symbol}
870       """
871       def __init__(self,arg,index):
872          """
873          initialization of wherePositive L{Symbol} with argument arg
874          @param arg: argument
875          @type arg: L{Symbol}.
876          @param index: defines index
877          @type index: C{slice} or C{int} or a C{tuple} of them
878          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
879          @raises ValueError: if a step is given
880          """
881          if not isinstance(index,tuple): index=(index,)
882          if len(index)>arg.getRank():
883               raise IndexError,"GetSlice_Symbol: index out of range."
884          sh=()
885          index2=()
886          for i in range(len(index)):
887             ix=index[i]
888             if isinstance(ix,int):
889                if ix<0 or ix>=arg.getShape()[i]:
890                   raise ValueError,"GetSlice_Symbol: index out of range."
891                index2=index2+(ix,)
892             else:
893               if not ix.step==None:
894                 raise ValueError,"GetSlice_Symbol: steping is not supported."
895               if ix.start==None:
896                  s=0
897               else:
898                  s=ix.start
899               if ix.stop==None:
900                  e=arg.getShape()[i]
901               else:
902                  e=ix.stop
903                  if e>arg.getShape()[i]:
904                     raise IndexError,"GetSlice_Symbol: index out of range."
905               index2=index2+(slice(s,e),)
906               if e>s:
907                   sh=sh+(e-s,)
908               elif s>e:
909                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
910          for i in range(len(index),arg.getRank()):
911              index2=index2+(slice(0,arg.getShape()[i]),)
912              sh=sh+(arg.getShape()[i],)
913          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
914    
915       def getMyCode(self,argstrs,format="escript"):
916          """
917          returns a program code that can be used to evaluate the symbol.
918    
919          @param argstrs: gives for each argument a string representing the argument for the evaluation.
920          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
921          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
922          @type format: C{str}
923          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
924          @rtype: C{str}
925          @raise NotImplementedError: if the requested format is not available
926          """
927          if format=="escript" or format=="str"  or format=="text":
928             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
929          else:
930             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
931    
932       def substitute(self,argvals):
933          """
934          assigns new values to symbols in the definition of the symbol.
935          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
936    
937          @param argvals: new values assigned to symbols
938          @type argvals: C{dict} with keywords of type L{Symbol}.
939          @return: result of the substitution process. Operations are executed as much as possible.
940          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
941          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
942          """
943          if argvals.has_key(self):
944             arg=argvals[self]
945             if self.isAppropriateValue(arg):
946                return arg
947             else:
948                raise TypeError,"%s: new value is not appropriate."%str(self)
949          else:
950             args=self.getSubstitutedArguments(argvals)
951             arg=args[0]
952             index=args[1]
953             return arg.__getitem__(index)
954    
955  def log10(arg):  def log10(arg):
956     """     """
957     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
958    
959     @param arg: argument     @param arg: argument
960     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
961     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
962     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
963     """     """
964     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 903  def wherePositive(arg): Line 980  def wherePositive(arg):
980    
981     @param arg: argument     @param arg: argument
982     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
983     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
984     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
985     """     """
986     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
987        if arg.rank==0:        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
988           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
989             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))  
990     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
991        return arg._wherePositive()        return arg._wherePositive()
992     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 953  class WherePositive_Symbol(DependendSymb Line 1026  class WherePositive_Symbol(DependendSymb
1026        @type format: C{str}        @type format: C{str}
1027        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1028        @rtype: C{str}        @rtype: C{str}
1029        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1030        """        """
1031        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1032            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 989  def whereNegative(arg): Line 1062  def whereNegative(arg):
1062    
1063     @param arg: argument     @param arg: argument
1064     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1065     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1066     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1067     """     """
1068     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1069        if arg.rank==0:        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1070           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1071             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))  
1072     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1073        return arg._whereNegative()        return arg._whereNegative()
1074     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1039  class WhereNegative_Symbol(DependendSymb Line 1108  class WhereNegative_Symbol(DependendSymb
1108        @type format: C{str}        @type format: C{str}
1109        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1110        @rtype: C{str}        @rtype: C{str}
1111        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1112        """        """
1113        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1114            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1075  def whereNonNegative(arg): Line 1144  def whereNonNegative(arg):
1144    
1145     @param arg: argument     @param arg: argument
1146     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1147     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1148     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1149     """     """
1150     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1151        if arg.rank==0:        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1152           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1153             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))  
1154     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1155        return arg._whereNonNegative()        return arg._whereNonNegative()
1156     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1109  def whereNonPositive(arg): Line 1174  def whereNonPositive(arg):
1174    
1175     @param arg: argument     @param arg: argument
1176     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1177     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1178     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1179     """     """
1180     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1181        if arg.rank==0:        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1182           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1183             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.  
1184     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1185        return arg._whereNonPositive()        return arg._whereNonPositive()
1186     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1145  def whereZero(arg,tol=0.): Line 1206  def whereZero(arg,tol=0.):
1206     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1207     @param tol: tolerance. values with absolute value less then tol are accepted as zero.     @param tol: tolerance. values with absolute value less then tol are accepted as zero.
1208     @type tol: C{float}     @type tol: C{float}
1209     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1210     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1211     """     """
1212     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1213        if arg.rank==0:        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1214           if abs(arg)<=tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1215             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
1216     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1217        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1218     elif isinstance(arg,float):     elif isinstance(arg,float):
1219        if abs(arg)<=tol:        if abs(arg)<=tol:
1220          return 1.          return 1.
# Line 1198  class WhereZero_Symbol(DependendSymbol): Line 1252  class WhereZero_Symbol(DependendSymbol):
1252        @type format: C{str}        @type format: C{str}
1253        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1254        @rtype: C{str}        @rtype: C{str}
1255        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1256        """        """
1257        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1258           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1232  def whereNonZero(arg,tol=0.): Line 1286  def whereNonZero(arg,tol=0.):
1286    
1287     @param arg: argument     @param arg: argument
1288     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1289     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1290     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1291     """     """
1292     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1293        if arg.rank==0:        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1294          if abs(arg)>tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1295             return numarray.array(1.)        return out
         else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
1296     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1297        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1298     elif isinstance(arg,float):     elif isinstance(arg,float):
1299        if abs(arg)>tol:        if abs(arg)>tol:
1300          return 1.          return 1.
# Line 1269  def sin(arg): Line 1316  def sin(arg):
1316    
1317     @param arg: argument     @param arg: argument
1318     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1319     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1320     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1321     """     """
1322     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1307  class Sin_Symbol(DependendSymbol): Line 1354  class Sin_Symbol(DependendSymbol):
1354        @type format: C{str}        @type format: C{str}
1355        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1356        @rtype: C{str}        @rtype: C{str}
1357        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1358        """        """
1359        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1360            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1359  def cos(arg): Line 1406  def cos(arg):
1406    
1407     @param arg: argument     @param arg: argument
1408     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1409     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1410     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1411     """     """
1412     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1397  class Cos_Symbol(DependendSymbol): Line 1444  class Cos_Symbol(DependendSymbol):
1444        @type format: C{str}        @type format: C{str}
1445        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1446        @rtype: C{str}        @rtype: C{str}
1447        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1448        """        """
1449        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1450            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1449  def tan(arg): Line 1496  def tan(arg):
1496    
1497     @param arg: argument     @param arg: argument
1498     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1499     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1500     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1501     """     """
1502     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1487  class Tan_Symbol(DependendSymbol): Line 1534  class Tan_Symbol(DependendSymbol):
1534        @type format: C{str}        @type format: C{str}
1535        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1536        @rtype: C{str}        @rtype: C{str}
1537        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1538        """        """
1539        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1540            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1539  def asin(arg): Line 1586  def asin(arg):
1586    
1587     @param arg: argument     @param arg: argument
1588     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1589     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1590     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1591     """     """
1592     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1577  class Asin_Symbol(DependendSymbol): Line 1624  class Asin_Symbol(DependendSymbol):
1624        @type format: C{str}        @type format: C{str}
1625        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1626        @rtype: C{str}        @rtype: C{str}
1627        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1628        """        """
1629        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1630            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1629  def acos(arg): Line 1676  def acos(arg):
1676    
1677     @param arg: argument     @param arg: argument
1678     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1679     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1680     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1681     """     """
1682     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1667  class Acos_Symbol(DependendSymbol): Line 1714  class Acos_Symbol(DependendSymbol):
1714        @type format: C{str}        @type format: C{str}
1715        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1716        @rtype: C{str}        @rtype: C{str}
1717        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1718        """        """
1719        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1720            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1719  def atan(arg): Line 1766  def atan(arg):
1766    
1767     @param arg: argument     @param arg: argument
1768     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1769     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1770     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1771     """     """
1772     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1757  class Atan_Symbol(DependendSymbol): Line 1804  class Atan_Symbol(DependendSymbol):
1804        @type format: C{str}        @type format: C{str}
1805        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1806        @rtype: C{str}        @rtype: C{str}
1807        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1808        """        """
1809        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1810            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1809  def sinh(arg): Line 1856  def sinh(arg):
1856    
1857     @param arg: argument     @param arg: argument
1858     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1859     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1860     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1861     """     """
1862     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1847  class Sinh_Symbol(DependendSymbol): Line 1894  class Sinh_Symbol(DependendSymbol):
1894        @type format: C{str}        @type format: C{str}
1895        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1896        @rtype: C{str}        @rtype: C{str}
1897        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1898        """        """
1899        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1900            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1899  def cosh(arg): Line 1946  def cosh(arg):
1946    
1947     @param arg: argument     @param arg: argument
1948     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1949     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1950     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1951     """     """
1952     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1937  class Cosh_Symbol(DependendSymbol): Line 1984  class Cosh_Symbol(DependendSymbol):
1984        @type format: C{str}        @type format: C{str}
1985        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1986        @rtype: C{str}        @rtype: C{str}
1987        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1988        """        """
1989        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1990            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1989  def tanh(arg): Line 2036  def tanh(arg):
2036    
2037     @param arg: argument     @param arg: argument
2038     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2039     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2040     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2041     """     """
2042     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2027  class Tanh_Symbol(DependendSymbol): Line 2074  class Tanh_Symbol(DependendSymbol):
2074        @type format: C{str}        @type format: C{str}
2075        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2076        @rtype: C{str}        @rtype: C{str}
2077        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2078        """        """
2079        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2080            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2079  def asinh(arg): Line 2126  def asinh(arg):
2126    
2127     @param arg: argument     @param arg: argument
2128     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2129     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2130     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2131     """     """
2132     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2117  class Asinh_Symbol(DependendSymbol): Line 2164  class Asinh_Symbol(DependendSymbol):
2164        @type format: C{str}        @type format: C{str}
2165        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2166        @rtype: C{str}        @rtype: C{str}
2167        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2168        """        """
2169        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2170            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2169  def acosh(arg): Line 2216  def acosh(arg):
2216    
2217     @param arg: argument     @param arg: argument
2218     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2219     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2220     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2221     """     """
2222     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2207  class Acosh_Symbol(DependendSymbol): Line 2254  class Acosh_Symbol(DependendSymbol):
2254        @type format: C{str}        @type format: C{str}
2255        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2256        @rtype: C{str}        @rtype: C{str}
2257        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2258        """        """
2259        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2260            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2259  def atanh(arg): Line 2306  def atanh(arg):
2306    
2307     @param arg: argument     @param arg: argument
2308     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2309     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2310     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2311     """     """
2312     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2297  class Atanh_Symbol(DependendSymbol): Line 2344  class Atanh_Symbol(DependendSymbol):
2344        @type format: C{str}        @type format: C{str}
2345        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2346        @rtype: C{str}        @rtype: C{str}
2347        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2348        """        """
2349        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2350            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2349  def exp(arg): Line 2396  def exp(arg):
2396    
2397     @param arg: argument     @param arg: argument
2398     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2399     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2400     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2401     """     """
2402     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2387  class Exp_Symbol(DependendSymbol): Line 2434  class Exp_Symbol(DependendSymbol):
2434        @type format: C{str}        @type format: C{str}
2435        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2436        @rtype: C{str}        @rtype: C{str}
2437        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2438        """        """
2439        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2440            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2439  def sqrt(arg): Line 2486  def sqrt(arg):
2486    
2487     @param arg: argument     @param arg: argument
2488     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2489     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2490     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2491     """     """
2492     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2477  class Sqrt_Symbol(DependendSymbol): Line 2524  class Sqrt_Symbol(DependendSymbol):
2524        @type format: C{str}        @type format: C{str}
2525        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2526        @rtype: C{str}        @rtype: C{str}
2527        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2528        """        """
2529        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2530            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2529  def log(arg): Line 2576  def log(arg):
2576    
2577     @param arg: argument     @param arg: argument
2578     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2579     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2580     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2581     """     """
2582     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2567  class Log_Symbol(DependendSymbol): Line 2614  class Log_Symbol(DependendSymbol):
2614        @type format: C{str}        @type format: C{str}
2615        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2616        @rtype: C{str}        @rtype: C{str}
2617        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2618        """        """
2619        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2620            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2619  def sign(arg): Line 2666  def sign(arg):
2666    
2667     @param arg: argument     @param arg: argument
2668     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2669     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2670     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2671     """     """
2672     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2667  class Abs_Symbol(DependendSymbol): Line 2714  class Abs_Symbol(DependendSymbol):
2714        @type format: C{str}        @type format: C{str}
2715        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2716        @rtype: C{str}        @rtype: C{str}
2717        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2718        """        """
2719        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2720            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2719  def minval(arg): Line 2766  def minval(arg):
2766    
2767     @param arg: argument     @param arg: argument
2768     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2769     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2770     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2771     """     """
2772     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2760  class Minval_Symbol(DependendSymbol): Line 2807  class Minval_Symbol(DependendSymbol):
2807        @type format: C{str}        @type format: C{str}
2808        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2809        @rtype: C{str}        @rtype: C{str}
2810        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2811        """        """
2812        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2813            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2796  def maxval(arg): Line 2843  def maxval(arg):
2843    
2844     @param arg: argument     @param arg: argument
2845     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2846     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2847     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2848     """     """
2849     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2837  class Maxval_Symbol(DependendSymbol): Line 2884  class Maxval_Symbol(DependendSymbol):
2884        @type format: C{str}        @type format: C{str}
2885        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2886        @rtype: C{str}        @rtype: C{str}
2887        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2888        """        """
2889        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2890            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2873  def length(arg): Line 2920  def length(arg):
2920    
2921     @param arg: argument     @param arg: argument
2922     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2923     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2924     """     """
2925     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2926    
2927    def trace(arg,axis_offset=0):
2928       """
2929       returns the trace of arg which the sum of arg[k,k] over k.
2930    
2931       @param arg: argument
2932       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2933       @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
2934                      C{axis_offset} and axis_offset+1 must be equal.
2935       @type axis_offset: C{int}
2936       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2937       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2938       """
2939       if isinstance(arg,numarray.NumArray):
2940          sh=arg.shape
2941          if len(sh)<2:
2942            raise ValueError,"rank of argument must be greater than 1"
2943          if axis_offset<0 or axis_offset>len(sh)-2:
2944            raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
2945          s1=1
2946          for i in range(axis_offset): s1*=sh[i]
2947          s2=1
2948          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2949          if not sh[axis_offset] == sh[axis_offset+1]:
2950            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2951          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
2952          out=numarray.zeros([s1,s2],numarray.Float64)
2953          for i1 in range(s1):
2954            for i2 in range(s2):
2955                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
2956          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2957          return out
2958       elif isinstance(arg,escript.Data):
2959          if arg.getRank()<2:
2960            raise ValueError,"rank of argument must be greater than 1"
2961          if axis_offset<0 or axis_offset>arg.getRank()-2:
2962            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
2963          s=list(arg.getShape())        
2964          if not s[axis_offset] == s[axis_offset+1]:
2965            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2966          return arg._trace(axis_offset)
2967       elif isinstance(arg,float):
2968          raise TypeError,"illegal argument type float."
2969       elif isinstance(arg,int):
2970          raise TypeError,"illegal argument type int."
2971       elif isinstance(arg,Symbol):
2972          return Trace_Symbol(arg,axis_offset)
2973       else:
2974          raise TypeError,"Unknown argument type."
2975    
2976    class Trace_Symbol(DependendSymbol):
2977       """
2978       L{Symbol} representing the result of the trace function
2979       """
2980       def __init__(self,arg,axis_offset=0):
2981          """
2982          initialization of trace L{Symbol} with argument arg
2983          @param arg: argument of function
2984          @type arg: L{Symbol}.
2985          @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
2986                      C{axis_offset} and axis_offset+1 must be equal.
2987          @type axis_offset: C{int}
2988          """
2989          if arg.getRank()<2:
2990            raise ValueError,"rank of argument must be greater than 1"
2991          if axis_offset<0 or axis_offset>arg.getRank()-2:
2992            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
2993          s=list(arg.getShape())        
2994          if not s[axis_offset] == s[axis_offset+1]:
2995            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2996          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
2997    
2998       def getMyCode(self,argstrs,format="escript"):
2999          """
3000          returns a program code that can be used to evaluate the symbol.
3001    
3002          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3003          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3004          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3005          @type format: C{str}
3006          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3007          @rtype: C{str}
3008          @raise NotImplementedError: if the requested format is not available
3009          """
3010          if format=="escript" or format=="str"  or format=="text":
3011             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3012          else:
3013             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3014    
3015       def substitute(self,argvals):
3016          """
3017          assigns new values to symbols in the definition of the symbol.
3018          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3019    
3020          @param argvals: new values assigned to symbols
3021          @type argvals: C{dict} with keywords of type L{Symbol}.
3022          @return: result of the substitution process. Operations are executed as much as possible.
3023          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3024          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3025          """
3026          if argvals.has_key(self):
3027             arg=argvals[self]
3028             if self.isAppropriateValue(arg):
3029                return arg
3030             else:
3031                raise TypeError,"%s: new value is not appropriate."%str(self)
3032          else:
3033             arg=self.getSubstitutedArguments(argvals)
3034             return trace(arg[0],axis_offset=arg[1])
3035    
3036       def diff(self,arg):
3037          """
3038          differential of this object
3039    
3040          @param arg: the derivative is calculated with respect to arg
3041          @type arg: L{escript.Symbol}
3042          @return: derivative with respect to C{arg}
3043          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3044          """
3045          if arg==self:
3046             return identity(self.getShape())
3047          else:
3048             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3049    
3050    def transpose(arg,axis_offset=None):
3051       """
3052       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3053    
3054       @param arg: argument
3055       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3056       @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.
3057                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3058       @type axis_offset: C{int}
3059       @return: transpose of arg
3060       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3061       """
3062       if isinstance(arg,numarray.NumArray):
3063          if axis_offset==None: axis_offset=int(arg.rank/2)
3064          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3065       elif isinstance(arg,escript.Data):
3066          r=arg.getRank()
3067          if axis_offset==None: axis_offset=int(r/2)
3068          if axis_offset<0 or axis_offset>r:
3069            raise ValueError,"axis_offset must be between 0 and %s"%r
3070          return arg._transpose(axis_offset)
3071       elif isinstance(arg,float):
3072          if not ( axis_offset==0 or axis_offset==None):
3073            raise ValueError,"axis_offset must be 0 for float argument"
3074          return arg
3075       elif isinstance(arg,int):
3076          if not ( axis_offset==0 or axis_offset==None):
3077            raise ValueError,"axis_offset must be 0 for int argument"
3078          return float(arg)
3079       elif isinstance(arg,Symbol):
3080          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3081          return Transpose_Symbol(arg,axis_offset)
3082       else:
3083          raise TypeError,"Unknown argument type."
3084    
3085    class Transpose_Symbol(DependendSymbol):
3086       """
3087       L{Symbol} representing the result of the transpose function
3088       """
3089       def __init__(self,arg,axis_offset=None):
3090          """
3091          initialization of transpose L{Symbol} with argument arg
3092    
3093          @param arg: argument of function
3094          @type arg: L{Symbol}.
3095          @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.
3096                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3097          @type axis_offset: C{int}
3098          """
3099          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3100          if axis_offset<0 or axis_offset>arg.getRank():
3101            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3102          s=arg.getShape()
3103          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3104    
3105       def getMyCode(self,argstrs,format="escript"):
3106          """
3107          returns a program code that can be used to evaluate the symbol.
3108    
3109          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3110          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3111          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3112          @type format: C{str}
3113          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3114          @rtype: C{str}
3115          @raise NotImplementedError: if the requested format is not available
3116          """
3117          if format=="escript" or format=="str"  or format=="text":
3118             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3119          else:
3120             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3121    
3122       def substitute(self,argvals):
3123          """
3124          assigns new values to symbols in the definition of the symbol.
3125          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3126    
3127          @param argvals: new values assigned to symbols
3128          @type argvals: C{dict} with keywords of type L{Symbol}.
3129          @return: result of the substitution process. Operations are executed as much as possible.
3130          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3131          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3132          """
3133          if argvals.has_key(self):
3134             arg=argvals[self]
3135             if self.isAppropriateValue(arg):
3136                return arg
3137             else:
3138                raise TypeError,"%s: new value is not appropriate."%str(self)
3139          else:
3140             arg=self.getSubstitutedArguments(argvals)
3141             return transpose(arg[0],axis_offset=arg[1])
3142    
3143       def diff(self,arg):
3144          """
3145          differential of this object
3146    
3147          @param arg: the derivative is calculated with respect to arg
3148          @type arg: L{escript.Symbol}
3149          @return: derivative with respect to C{arg}
3150          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3151          """
3152          if arg==self:
3153             return identity(self.getShape())
3154          else:
3155             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3156    
3157    def swap_axes(arg,axis0=0,axis1=1):
3158       """
3159       returns the swap of arg by swaping the components axis0 and axis1
3160    
3161       @param arg: argument
3162       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3163       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3164       @type axis0: C{int}
3165       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3166       @type axis1: C{int}
3167       @return: C{arg} with swaped components
3168       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3169       """
3170       if axis0 > axis1:
3171          axis0,axis1=axis1,axis0
3172       if isinstance(arg,numarray.NumArray):
3173          return numarray.swapaxes(arg,axis0,axis1)
3174       elif isinstance(arg,escript.Data):
3175          return arg._swap_axes(axis0,axis1)
3176       elif isinstance(arg,float):
3177          raise TyepError,"float argument is not supported."
3178       elif isinstance(arg,int):
3179          raise TyepError,"int argument is not supported."
3180       elif isinstance(arg,Symbol):
3181          return SwapAxes_Symbol(arg,axis0,axis1)
3182       else:
3183          raise TypeError,"Unknown argument type."
3184    
3185    class SwapAxes_Symbol(DependendSymbol):
3186       """
3187       L{Symbol} representing the result of the swap function
3188       """
3189       def __init__(self,arg,axis0=0,axis1=1):
3190          """
3191          initialization of swap L{Symbol} with argument arg
3192    
3193          @param arg: argument
3194          @type arg: L{Symbol}.
3195          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3196          @type axis0: C{int}
3197          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3198          @type axis1: C{int}
3199          """
3200          if arg.getRank()<2:
3201             raise ValueError,"argument must have at least rank 2."
3202          if axis0<0 or axis0>arg.getRank()-1:
3203             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3204          if axis1<0 or axis1>arg.getRank()-1:
3205             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3206          if axis0 == axis1:
3207             raise ValueError,"axis indices must be different."
3208          if axis0 > axis1:
3209             axis0,axis1=axis1,axis0
3210          s=arg.getShape()
3211          s_out=[]
3212          for i in range(len(s)):
3213             if i == axis0:
3214                s_out.append(s[axis1])
3215             elif i == axis1:
3216                s_out.append(s[axis0])
3217             else:
3218                s_out.append(s[i])
3219          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3220    
3221       def getMyCode(self,argstrs,format="escript"):
3222          """
3223          returns a program code that can be used to evaluate the symbol.
3224    
3225          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3226          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3227          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3228          @type format: C{str}
3229          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3230          @rtype: C{str}
3231          @raise NotImplementedError: if the requested format is not available
3232          """
3233          if format=="escript" or format=="str"  or format=="text":
3234             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3235          else:
3236             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3237    
3238       def substitute(self,argvals):
3239          """
3240          assigns new values to symbols in the definition of the symbol.
3241          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3242    
3243          @param argvals: new values assigned to symbols
3244          @type argvals: C{dict} with keywords of type L{Symbol}.
3245          @return: result of the substitution process. Operations are executed as much as possible.
3246          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3247          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3248          """
3249          if argvals.has_key(self):
3250             arg=argvals[self]
3251             if self.isAppropriateValue(arg):
3252                return arg
3253             else:
3254                raise TypeError,"%s: new value is not appropriate."%str(self)
3255          else:
3256             arg=self.getSubstitutedArguments(argvals)
3257             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3258    
3259       def diff(self,arg):
3260          """
3261          differential of this object
3262    
3263          @param arg: the derivative is calculated with respect to arg
3264          @type arg: L{escript.Symbol}
3265          @return: derivative with respect to C{arg}
3266          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3267          """
3268          if arg==self:
3269             return identity(self.getShape())
3270          else:
3271             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3272    
3273    def symmetric(arg):
3274        """
3275        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3276    
3277        @param arg: square matrix. Must have rank 2 or 4 and be square.
3278        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3279        @return: symmetric part of arg
3280        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3281        """
3282        if isinstance(arg,numarray.NumArray):
3283          if arg.rank==2:
3284            if not (arg.shape[0]==arg.shape[1]):
3285               raise ValueError,"argument must be square."
3286          elif arg.rank==4:
3287            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3288               raise ValueError,"argument must be square."
3289          else:
3290            raise ValueError,"rank 2 or 4 is required."
3291          return (arg+transpose(arg))/2
3292        elif isinstance(arg,escript.Data):
3293          if arg.getRank()==2:
3294            if not (arg.getShape()[0]==arg.getShape()[1]):
3295               raise ValueError,"argument must be square."
3296            return arg._symmetric()
3297          elif arg.getRank()==4:
3298            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3299               raise ValueError,"argument must be square."
3300            return arg._symmetric()
3301          else:
3302            raise ValueError,"rank 2 or 4 is required."
3303        elif isinstance(arg,float):
3304          return arg
3305        elif isinstance(arg,int):
3306          return float(arg)
3307        elif isinstance(arg,Symbol):
3308          if arg.getRank()==2:
3309            if not (arg.getShape()[0]==arg.getShape()[1]):
3310               raise ValueError,"argument must be square."
3311          elif arg.getRank()==4:
3312            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3313               raise ValueError,"argument must be square."
3314          else:
3315            raise ValueError,"rank 2 or 4 is required."
3316          return (arg+transpose(arg))/2
3317        else:
3318          raise TypeError,"symmetric: Unknown argument type."
3319    
3320    def nonsymmetric(arg):
3321        """
3322        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3323    
3324        @param arg: square matrix. Must have rank 2 or 4 and be square.
3325        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3326        @return: nonsymmetric part of arg
3327        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3328        """
3329        if isinstance(arg,numarray.NumArray):
3330          if arg.rank==2:
3331            if not (arg.shape[0]==arg.shape[1]):
3332               raise ValueError,"nonsymmetric: argument must be square."
3333          elif arg.rank==4:
3334            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3335               raise ValueError,"nonsymmetric: argument must be square."
3336          else:
3337            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3338          return (arg-transpose(arg))/2
3339        elif isinstance(arg,escript.Data):
3340          if arg.getRank()==2:
3341            if not (arg.getShape()[0]==arg.getShape()[1]):
3342               raise ValueError,"argument must be square."
3343            return arg._nonsymmetric()
3344          elif arg.getRank()==4:
3345            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3346               raise ValueError,"argument must be square."
3347            return arg._nonsymmetric()
3348          else:
3349            raise ValueError,"rank 2 or 4 is required."
3350        elif isinstance(arg,float):
3351          return arg
3352        elif isinstance(arg,int):
3353          return float(arg)
3354        elif isinstance(arg,Symbol):
3355          if arg.getRank()==2:
3356            if not (arg.getShape()[0]==arg.getShape()[1]):
3357               raise ValueError,"nonsymmetric: argument must be square."
3358          elif arg.getRank()==4:
3359            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3360               raise ValueError,"nonsymmetric: argument must be square."
3361          else:
3362            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3363          return (arg-transpose(arg))/2
3364        else:
3365          raise TypeError,"nonsymmetric: Unknown argument type."
3366    
3367    def inverse(arg):
3368        """
3369        returns the inverse of the square matrix arg.
3370    
3371        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3372        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3373        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3374        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3375        @note: for L{escript.Data} objects the dimension is restricted to 3.
3376        """
3377        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3378        if isinstance(arg,numarray.NumArray):
3379          return numarray.linear_algebra.inverse(arg)
3380        elif isinstance(arg,escript.Data):
3381          return escript_inverse(arg)
3382        elif isinstance(arg,float):
3383          return 1./arg
3384        elif isinstance(arg,int):
3385          return 1./float(arg)
3386        elif isinstance(arg,Symbol):
3387          return Inverse_Symbol(arg)
3388        else:
3389          raise TypeError,"inverse: Unknown argument type."
3390    
3391    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3392          "arg is a Data objects!!!"
3393          if not arg.getRank()==2:
3394            raise ValueError,"escript_inverse: argument must have rank 2"
3395          s=arg.getShape()      
3396          if not s[0] == s[1]:
3397            raise ValueError,"escript_inverse: argument must be a square matrix."
3398          out=escript.Data(0.,s,arg.getFunctionSpace())
3399          if s[0]==1:
3400              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3401                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3402              out[0,0]=1./arg[0,0]
3403          elif s[0]==2:
3404              A11=arg[0,0]
3405              A12=arg[0,1]
3406              A21=arg[1,0]
3407              A22=arg[1,1]
3408              D = A11*A22-A12*A21
3409              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3410                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3411              D=1./D
3412              out[0,0]= A22*D
3413              out[1,0]=-A21*D
3414              out[0,1]=-A12*D
3415              out[1,1]= A11*D
3416          elif s[0]==3:
3417              A11=arg[0,0]
3418              A21=arg[1,0]
3419              A31=arg[2,0]
3420              A12=arg[0,1]
3421              A22=arg[1,1]
3422              A32=arg[2,1]
3423              A13=arg[0,2]
3424              A23=arg[1,2]
3425              A33=arg[2,2]
3426              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3427              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3428                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3429              D=1./D
3430              out[0,0]=(A22*A33-A23*A32)*D
3431              out[1,0]=(A31*A23-A21*A33)*D
3432              out[2,0]=(A21*A32-A31*A22)*D
3433              out[0,1]=(A13*A32-A12*A33)*D
3434              out[1,1]=(A11*A33-A31*A13)*D
3435              out[2,1]=(A12*A31-A11*A32)*D
3436              out[0,2]=(A12*A23-A13*A22)*D
3437              out[1,2]=(A13*A21-A11*A23)*D
3438              out[2,2]=(A11*A22-A12*A21)*D
3439          else:
3440             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3441          return out
3442    
3443    class Inverse_Symbol(DependendSymbol):
3444       """
3445       L{Symbol} representing the result of the inverse function
3446       """
3447       def __init__(self,arg):
3448          """
3449          initialization of inverse L{Symbol} with argument arg
3450          @param arg: argument of function
3451          @type arg: L{Symbol}.
3452          """
3453          if not arg.getRank()==2:
3454            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3455          s=arg.getShape()
3456          if not s[0] == s[1]:
3457            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3458          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3459    
3460       def getMyCode(self,argstrs,format="escript"):
3461          """
3462          returns a program code that can be used to evaluate the symbol.
3463    
3464          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3465          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3466          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3467          @type format: C{str}
3468          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3469          @rtype: C{str}
3470          @raise NotImplementedError: if the requested format is not available
3471          """
3472          if format=="escript" or format=="str"  or format=="text":
3473             return "inverse(%s)"%argstrs[0]
3474          else:
3475             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3476    
3477       def substitute(self,argvals):
3478          """
3479          assigns new values to symbols in the definition of the symbol.
3480          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3481    
3482          @param argvals: new values assigned to symbols
3483          @type argvals: C{dict} with keywords of type L{Symbol}.
3484          @return: result of the substitution process. Operations are executed as much as possible.
3485          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3486          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3487          """
3488          if argvals.has_key(self):
3489             arg=argvals[self]
3490             if self.isAppropriateValue(arg):
3491                return arg
3492             else:
3493                raise TypeError,"%s: new value is not appropriate."%str(self)
3494          else:
3495             arg=self.getSubstitutedArguments(argvals)
3496             return inverse(arg[0])
3497    
3498       def diff(self,arg):
3499          """
3500          differential of this object
3501    
3502          @param arg: the derivative is calculated with respect to arg
3503          @type arg: L{escript.Symbol}
3504          @return: derivative with respect to C{arg}
3505          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3506          """
3507          if arg==self:
3508             return identity(self.getShape())
3509          else:
3510             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3511    
3512    def eigenvalues(arg):
3513        """
3514        returns the eigenvalues of the square matrix arg.
3515    
3516        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3517                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3518        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3519        @return: the eigenvalues in increasing order.
3520        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3521        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3522        """
3523        if isinstance(arg,numarray.NumArray):
3524          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3525          out.sort()
3526          return out
3527        elif isinstance(arg,escript.Data):
3528          return arg._eigenvalues()
3529        elif isinstance(arg,Symbol):
3530          if not arg.getRank()==2:
3531            raise ValueError,"eigenvalues: argument must have rank 2"
3532          s=arg.getShape()      
3533          if not s[0] == s[1]:
3534            raise ValueError,"eigenvalues: argument must be a square matrix."
3535          if s[0]==1:
3536              return arg[0]
3537          elif s[0]==2:
3538              arg1=symmetric(arg)
3539              A11=arg1[0,0]
3540              A12=arg1[0,1]
3541              A22=arg1[1,1]
3542              trA=(A11+A22)/2.
3543              A11-=trA
3544              A22-=trA
3545              s=sqrt(A12**2-A11*A22)
3546              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3547          elif s[0]==3:
3548              arg1=symmetric(arg)
3549              A11=arg1[0,0]
3550              A12=arg1[0,1]
3551              A22=arg1[1,1]
3552              A13=arg1[0,2]
3553              A23=arg1[1,2]
3554              A33=arg1[2,2]
3555              trA=(A11+A22+A33)/3.
3556              A11-=trA
3557              A22-=trA
3558              A33-=trA
3559              A13_2=A13**2
3560              A23_2=A23**2
3561              A12_2=A12**2
3562              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3563              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3564              sq_p=sqrt(p/3.)
3565              alpha_3=acos(clip(-q*(sq_p+whereZero(p,0.)*1.e-15)**(-3.)/2.,-1.,1.))/3.  # whereZero is protection against divison by zero
3566              sq_p*=2.
3567              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3568               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3569               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3570              return trA+sq_p*f
3571          else:
3572             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3573        elif isinstance(arg,float):
3574          return arg
3575        elif isinstance(arg,int):
3576          return float(arg)
3577        else:
3578          raise TypeError,"eigenvalues: Unknown argument type."
3579    
3580    def eigenvalues_and_eigenvectors(arg):
3581        """
3582        returns the eigenvalues and eigenvectors of the square matrix arg.
3583    
3584        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3585                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3586        @type arg: L{escript.Data}
3587        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3588                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3589                 the eigenvector coresponding to the i-th eigenvalue.
3590        @rtype: L{tuple} of L{escript.Data}.
3591        @note: The dimension is restricted to 3.
3592        """
3593        if isinstance(arg,numarray.NumArray):
3594          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3595        elif isinstance(arg,escript.Data):
3596          return arg._eigenvalues_and_eigenvectors()
3597        elif isinstance(arg,Symbol):
3598          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3599        elif isinstance(arg,float):
3600          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3601        elif isinstance(arg,int):
3602          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3603        else:
3604          raise TypeError,"eigenvalues: Unknown argument type."
3605  #=======================================================  #=======================================================
3606  #  Binary operations:  #  Binary operations:
3607  #=======================================================  #=======================================================
# Line 2936  class Add_Symbol(DependendSymbol): Line 3661  class Add_Symbol(DependendSymbol):
3661        @type format: C{str}        @type format: C{str}
3662        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3663        @rtype: C{str}        @rtype: C{str}
3664        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3665        """        """
3666        if format=="str" or format=="text":        if format=="str" or format=="text":
3667           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 2995  def mult(arg0,arg1): Line 3720  def mult(arg0,arg1):
3720         """         """
3721         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3722         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3723            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3724         else:         else:
3725            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3726                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3035  class Mult_Symbol(DependendSymbol): Line 3760  class Mult_Symbol(DependendSymbol):
3760        @type format: C{str}        @type format: C{str}
3761        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3762        @rtype: C{str}        @rtype: C{str}
3763        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3764        """        """
3765        if format=="str" or format=="text":        if format=="str" or format=="text":
3766           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3095  def quotient(arg0,arg1): Line 3820  def quotient(arg0,arg1):
3820         """         """
3821         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3822         if testForZero(args[0]):         if testForZero(args[0]):
3823            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3824         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3825            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3826               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3140  class Quotient_Symbol(DependendSymbol): Line 3865  class Quotient_Symbol(DependendSymbol):
3865        @type format: C{str}        @type format: C{str}
3866        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3867        @rtype: C{str}        @rtype: C{str}
3868        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3869        """        """
3870        if format=="str" or format=="text":        if format=="str" or format=="text":
3871           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3201  def power(arg0,arg1): Line 3926  def power(arg0,arg1):
3926         """         """
3927         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3928         if testForZero(args[0]):         if testForZero(args[0]):
3929            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3930         elif testForZero(args[1]):         elif testForZero(args[1]):
3931            return numarray.ones(args[0],numarray.Float)            return numarray.ones(pokeShape(args[1]),numarray.Float64)
3932         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3933            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3934         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):
# Line 3244  class Power_Symbol(DependendSymbol): Line 3969  class Power_Symbol(DependendSymbol):
3969        @type format: C{str}        @type format: C{str}
3970        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3971        @rtype: C{str}        @rtype: C{str}
3972        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3973        """        """
3974        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
3975           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3304  def maximum(*args): Line 4029  def maximum(*args):
4029         if out==None:         if out==None:
4030            out=a            out=a
4031         else:         else:
4032            m=whereNegative(out-a)            diff=add(a,-out)
4033            out=m*a+(1.-m)*out            out=add(out,mult(wherePositive(diff),diff))
4034      return out      return out
4035        
4036  def minimum(*arg):  def minimum(*args):
4037      """      """
4038      the minimum over arguments args      the minimum over arguments args
4039    
# Line 3322  def minimum(*arg): Line 4047  def minimum(*arg):
4047         if out==None:         if out==None:
4048            out=a            out=a
4049         else:         else:
4050            m=whereNegative(out-a)            diff=add(a,-out)
4051            out=m*out+(1.-m)*a            out=add(out,mult(whereNegative(diff),diff))
4052      return out      return out
4053    
4054    def clip(arg,minval=0.,maxval=1.):
4055        """
4056        cuts the values of arg between minval and maxval
4057    
4058        @param arg: argument
4059        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4060        @param minval: lower range
4061        @type minval: C{float}
4062        @param maxval: upper range
4063        @type maxval: C{float}
4064        @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4065                 less then maxval are unchanged.
4066        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4067        @raise ValueError: if minval>maxval
4068        """
4069        if minval>maxval:
4070           raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4071        return minimum(maximum(minval,arg),maxval)
4072    
4073        
4074  def inner(arg0,arg1):  def inner(arg0,arg1):
4075      """      """
# Line 3340  def inner(arg0,arg1): Line 4085  def inner(arg0,arg1):
4085      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4086      @param arg1: second argument      @param arg1: second argument
4087      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4088      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4089      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input
4090      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4091      """      """
# Line 3348  def inner(arg0,arg1): Line 4093  def inner(arg0,arg1):
4093      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4094      if not sh0==sh1:      if not sh0==sh1:
4095          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4096      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4097    
4098    def outer(arg0,arg1):
4099        """
4100        the outer product of the two argument:
4101    
4102        out[t,s]=arg0[t]*arg1[s]
4103    
4104        where
4105    
4106            - s runs through arg0.Shape
4107            - t runs through arg1.Shape
4108    
4109        @param arg0: first argument
4110        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4111        @param arg1: second argument
4112        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4113        @return: the outer product of arg0 and arg1 at each data point
4114        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4115        """
4116        return generalTensorProduct(arg0,arg1,axis_offset=0)
4117    
4118  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4119      """      """
4120        see L{matrix_mult}
4121        """
4122        return matrix_mult(arg0,arg1)
4123    
4124    def matrix_mult(arg0,arg1):
4125        """
4126      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4127    
4128      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4129    
4130            or      or
4131    
4132      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4133    
4134      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.
4135    
4136      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4137      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3376  def matrixmult(arg0,arg1): Line 4147  def matrixmult(arg0,arg1):
4147          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4148      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4149          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4150      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4151    
4152  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4153      """      """
4154      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  
4155      """      """
4156      return generalTensorProduct(arg0,arg1,offset=0)      return tensor_mult(arg0,arg1)
4157    
4158    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4159      """      """
4160      the tensor product of the two argument:      the tensor product of the two argument:
   
4161            
4162      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4163    
4164      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4165    
4166                   or      or
4167    
4168      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4169    
# Line 3415  def tensormult(arg0,arg1): Line 4172  def tensormult(arg0,arg1):
4172    
4173      out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2,s3]      out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2,s3]
4174                                
4175                   or      or
4176    
4177      out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2]      out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2]
4178    
4179                   or      or
4180    
4181      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]
4182    
4183      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  
4184      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.
4185    
4186      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4187      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3436  def tensormult(arg0,arg1): Line 4193  def tensormult(arg0,arg1):
4193      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
4194      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4195      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4196         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4197      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4198         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4199      else:      else:
4200          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4201    
4202  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4203      """      """
4204      generalized tensor product      generalized tensor product
4205    
4206      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4207    
4208      where s runs through arg0.Shape[:arg0.Rank-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[offset:]  
4209    
4210      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]
4211      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4212            - t runs through arg1.Shape[axis_offset:]
4213    
4214      @param arg0: first argument      @param arg0: first argument
4215      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4216      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4217      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4218      @return: the general tensor product of arg0 and arg1 at each data point.      @return: the general tensor product of arg0 and arg1 at each data point.
4219      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
# Line 3467  def generalTensorProduct(arg0,arg1,offse Line 4223  def generalTensorProduct(arg0,arg1,offse
4223      # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols      # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4224      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4225         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4226             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4227         else:         else:
4228             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4229                 raise ValueError,"generalTensorProduct: dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)                 raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4230             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4231             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4232             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4233             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4234             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4235             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4236             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4237             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4238             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4239             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4240             for i0 in range(d0):             for i0 in range(d0):
4241                      for i1 in range(d1):                      for i1 in range(d1):
4242                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4243             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4244             return out             return out
4245      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4246         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4247             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4248         else:         else:
4249             return escript_generalTensorProduct(arg0,arg1,offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,offset)             return escript_generalTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4250      else:            else:      
4251         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4252                                    
4253  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4254     """     """
4255     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4256     """     """
4257     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4258         """         """
4259         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4260    
4261         @param arg0: numerator         @param arg0: first argument
4262         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4263         @param arg1: denominator         @param arg1: second argument
4264         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4265         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4266         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4267         """         """
4268         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
4269         sh_arg1=pokeShape(arg1)         sh_arg1=pokeShape(arg1)
4270         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4271         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4272         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4273         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4274         if not sh01==sh10:         if not sh01==sh10:
4275             raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)             raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4276         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,offset])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4277    
4278     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4279        """        """
# Line 3529  class GeneralTensorProduct_Symbol(Depend Line 4285  class GeneralTensorProduct_Symbol(Depend
4285        @type format: C{str}        @type format: C{str}
4286        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4287        @rtype: C{str}        @rtype: C{str}
4288        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4289        """        """
4290        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4291           return "generalTensorProduct(%s,%s,offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4292        else:        else:
4293           raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)           raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4294    
# Line 3557  class GeneralTensorProduct_Symbol(Depend Line 4313  class GeneralTensorProduct_Symbol(Depend
4313           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4314           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4315    
4316  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProductNew(arg0,arg1,axis_offset):
4317      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4318      # calculate the return shape:      # calculate the return shape:
4319      shape0=arg0.getShape()[:arg0.getRank()-offset]      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4320      shape01=arg0.getShape()[arg0.getRank()-offset:]      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4321      shape10=arg1.getShape()[:offset]      shape10=arg1.getShape()[:axis_offset]
4322      shape1=arg1.getShape()[offset:]      shape1=arg1.getShape()[axis_offset:]
4323      if not shape01==shape10:      if not shape01==shape10:
4324          raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)          raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4325        # Figure out which functionspace to use (look at where operator+ is defined maybe in BinaryOp.h to get the logic for this)
4326        # fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4327        out=GeneralTensorProduct(arg0, arg1, axis_offset)
4328        return out
4329    
4330    def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct
4331        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4332        # calculate the return shape:
4333        shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4334        shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4335        shape10=arg1.getShape()[:axis_offset]
4336        shape1=arg1.getShape()[axis_offset:]
4337        if not shape01==shape10:
4338            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)
4339    
4340        # whatr function space should be used? (this here is not good!)
4341        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4342      # create return value:      # create return value:
4343      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())      out=escript.Data(0.,tuple(shape0+shape1),fs)
4344      #      #
4345      s0=[[]]      s0=[[]]
4346      for k in shape0:      for k in shape0:
# Line 3591  def escript_generalTensorProduct(arg0,ar Line 4363  def escript_generalTensorProduct(arg0,ar
4363    
4364      for i0 in s0:      for i0 in s0:
4365         for i1 in s1:         for i1 in s1:
4366           s=escript.Scalar(0.,arg0.getFunctionSpace())           s=escript.Scalar(0.,fs)
4367           for i01 in s01:           for i01 in s01:
4368              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))
4369           out.__setitem__(tuple(i0+i1),s)           out.__setitem__(tuple(i0+i1),s)
4370      return out      return out
4371    
4372    
4373    def transposed_matrix_mult(arg0,arg1):
4374        """
4375        transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4376    
4377        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
4378    
4379        or
4380    
4381        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4382    
4383        The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4384    
4385        The first dimension of arg0 and arg1 must match.
4386    
4387        @param arg0: first argument of rank 2
4388        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4389        @param arg1: second argument of at least rank 1
4390        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4391        @return: the product of the transposed of arg0 and arg1 at each data point
4392        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4393        @raise ValueError: if the shapes of the arguments are not appropriate
4394        """
4395        sh0=pokeShape(arg0)
4396        sh1=pokeShape(arg1)
4397        if not len(sh0)==2 :
4398            raise ValueError,"first argument must have rank 2"
4399        if not len(sh1)==2 and not len(sh1)==1:
4400            raise ValueError,"second argument must have rank 1 or 2"
4401        return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4402    
4403    def transposed_tensor_mult(arg0,arg1):
4404        """
4405        the tensor product of the transposed of the first and the second argument
4406        
4407        for arg0 of rank 2 this is
4408    
4409        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4410    
4411        or
4412    
4413        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4414    
4415      
4416        and for arg0 of rank 4 this is
4417    
4418        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4419                  
4420        or
4421    
4422        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4423    
4424        or
4425    
4426        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4427    
4428        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4429        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4430    
4431        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4432    
4433        @param arg0: first argument of rank 2 or 4
4434        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4435        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4436        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4437        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4438        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4439        """
4440        sh0=pokeShape(arg0)
4441        sh1=pokeShape(arg1)
4442        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4443           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4444        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4445           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4446        else:
4447            raise ValueError,"first argument must have rank 2 or 4"
4448    
4449    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4450        """
4451        generalized tensor product of transposed of arg0 and arg1:
4452    
4453        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4454    
4455        where
4456    
4457            - s runs through arg0.Shape[axis_offset:]
4458            - r runs trough arg0.Shape[:axis_offset]
4459            - t runs through arg1.Shape[axis_offset:]
4460    
4461        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4462        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4463    
4464        @param arg0: first argument
4465        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4466        @param arg1: second argument
4467        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4468        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4469        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4470        """
4471        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4472        arg0,arg1=matchType(arg0,arg1)
4473        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4474        if isinstance(arg0,numarray.NumArray):
4475           if isinstance(arg1,Symbol):
4476               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4477           else:
4478               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4479                   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)
4480               arg0_c=arg0.copy()
4481               arg1_c=arg1.copy()
4482               sh0,sh1=arg0.shape,arg1.shape
4483               d0,d1,d01=1,1,1
4484               for i in sh0[axis_offset:]: d0*=i
4485               for i in sh1[axis_offset:]: d1*=i
4486               for i in sh0[:axis_offset]: d01*=i
4487               arg0_c.resize((d01,d0))
4488               arg1_c.resize((d01,d1))
4489               out=numarray.zeros((d0,d1),numarray.Float64)
4490               for i0 in range(d0):
4491                        for i1 in range(d1):
4492                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4493               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4494               return out
4495        elif isinstance(arg0,escript.Data):
4496           if isinstance(arg1,Symbol):
4497               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4498           else:
4499               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4500        else:      
4501           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4502                    
4503    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4504       """
4505       Symbol representing the general tensor product of the transposed of arg0 and arg1
4506       """
4507       def __init__(self,arg0,arg1,axis_offset=0):
4508           """
4509           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4510    
4511           @param arg0: first argument
4512           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4513           @param arg1: second argument
4514           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4515           @raise ValueError: inconsistent dimensions of arguments.
4516           @note: if both arguments have a spatial dimension, they must equal.
4517           """
4518           sh_arg0=pokeShape(arg0)
4519           sh_arg1=pokeShape(arg1)
4520           sh01=sh_arg0[:axis_offset]
4521           sh10=sh_arg1[:axis_offset]
4522           sh0=sh_arg0[axis_offset:]
4523           sh1=sh_arg1[axis_offset:]
4524           if not sh01==sh10:
4525               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)
4526           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4527    
4528       def getMyCode(self,argstrs,format="escript"):
4529          """
4530          returns a program code that can be used to evaluate the symbol.
4531    
4532          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4533          @type argstrs: C{list} of length 2 of C{str}.
4534          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4535          @type format: C{str}
4536          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4537          @rtype: C{str}
4538          @raise NotImplementedError: if the requested format is not available
4539          """
4540          if format=="escript" or format=="str" or format=="text":
4541             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4542          else:
4543             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4544    
4545       def substitute(self,argvals):
4546          """
4547          assigns new values to symbols in the definition of the symbol.
4548          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4549    
4550          @param argvals: new values assigned to symbols
4551          @type argvals: C{dict} with keywords of type L{Symbol}.
4552          @return: result of the substitution process. Operations are executed as much as possible.
4553          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4554          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4555          """
4556          if argvals.has_key(self):
4557             arg=argvals[self]
4558             if self.isAppropriateValue(arg):
4559                return arg
4560             else:
4561                raise TypeError,"%s: new value is not appropriate."%str(self)
4562          else:
4563             args=self.getSubstitutedArguments(argvals)
4564             return generalTransposedTensorProduct(args[0],args[1],args[2])
4565    
4566    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4567        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4568        # calculate the return shape:
4569        shape01=arg0.getShape()[:axis_offset]
4570        shape10=arg1.getShape()[:axis_offset]
4571        shape0=arg0.getShape()[axis_offset:]
4572        shape1=arg1.getShape()[axis_offset:]
4573        if not shape01==shape10:
4574            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)
4575    
4576        # whatr function space should be used? (this here is not good!)
4577        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4578        # create return value:
4579        out=escript.Data(0.,tuple(shape0+shape1),fs)
4580        #
4581        s0=[[]]
4582        for k in shape0:
4583              s=[]
4584              for j in s0:
4585                    for i in range(k): s.append(j+[slice(i,i)])
4586              s0=s
4587        s1=[[]]
4588        for k in shape1:
4589              s=[]
4590              for j in s1:
4591                    for i in range(k): s.append(j+[slice(i,i)])
4592              s1=s
4593        s01=[[]]
4594        for k in shape01:
4595              s=[]
4596              for j in s01:
4597                    for i in range(k): s.append(j+[slice(i,i)])
4598              s01=s
4599    
4600        for i0 in s0:
4601           for i1 in s1:
4602             s=escript.Scalar(0.,fs)
4603             for i01 in s01:
4604                s+=arg0.__getitem__(tuple(i01+i0))*arg1.__getitem__(tuple(i01+i1))
4605             out.__setitem__(tuple(i0+i1),s)
4606        return out
4607    
4608    
4609    def matrix_transposed_mult(arg0,arg1):
4610        """
4611        matrix-transposed(matrix) product of the two argument:
4612    
4613        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4614    
4615        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4616    
4617        The last dimensions of arg0 and arg1 must match.
4618    
4619        @param arg0: first argument of rank 2
4620        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4621        @param arg1: second argument of rank 2
4622        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4623        @return: the product of arg0 and the transposed of arg1 at each data point
4624        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4625        @raise ValueError: if the shapes of the arguments are not appropriate
4626        """
4627        sh0=pokeShape(arg0)
4628        sh1=pokeShape(arg1)
4629        if not len(sh0)==2 :
4630            raise ValueError,"first argument must have rank 2"
4631        if not len(sh1)==2 and not len(sh1)==1:
4632            raise ValueError,"second argument must have rank 1 or 2"
4633        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4634    
4635    def tensor_transposed_mult(arg0,arg1):
4636        """
4637        the tensor product of the first and the transpose of the second argument
4638        
4639        for arg0 of rank 2 this is
4640    
4641        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4642    
4643        and for arg0 of rank 4 this is
4644    
4645        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4646                  
4647        or
4648    
4649        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4650    
4651        In the first case the the second dimension of arg0 and arg1 must match and  
4652        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4653    
4654        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4655    
4656        @param arg0: first argument of rank 2 or 4
4657        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4658        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4659        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4660        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4661        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4662        """
4663        sh0=pokeShape(arg0)
4664        sh1=pokeShape(arg1)
4665        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4666           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4667        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4668           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4669        else:
4670            raise ValueError,"first argument must have rank 2 or 4"
4671    
4672    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4673        """
4674        generalized tensor product of transposed of arg0 and arg1:
4675    
4676        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4677    
4678        where
4679    
4680            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4681            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4682            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4683    
4684        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4685        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4686    
4687        @param arg0: first argument
4688        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4689        @param arg1: second argument
4690        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4691        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4692        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4693        """
4694        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4695        arg0,arg1=matchType(arg0,arg1)
4696        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4697        if isinstance(arg0,numarray.NumArray):
4698           if isinstance(arg1,Symbol):
4699               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4700           else:
4701               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4702                   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)
4703               arg0_c=arg0.copy()
4704               arg1_c=arg1.copy()
4705               sh0,sh1=arg0.shape,arg1.shape
4706               d0,d1,d01=1,1,1
4707               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4708               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4709               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4710               arg0_c.resize((d0,d01))
4711               arg1_c.resize((d1,d01))
4712               out=numarray.zeros((d0,d1),numarray.Float64)
4713               for i0 in range(d0):
4714                        for i1 in range(d1):
4715                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4716               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4717               return out
4718        elif isinstance(arg0,escript.Data):
4719           if isinstance(arg1,Symbol):
4720               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4721           else:
4722               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4723        else:      
4724           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4725                    
4726    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4727       """
4728       Symbol representing the general tensor product of arg0 and the transpose of arg1
4729       """
4730       def __init__(self,arg0,arg1,axis_offset=0):
4731           """
4732           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4733    
4734           @param arg0: first argument
4735           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4736           @param arg1: second argument
4737           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4738           @raise ValueError: inconsistent dimensions of arguments.
4739           @note: if both arguments have a spatial dimension, they must equal.
4740           """
4741           sh_arg0=pokeShape(arg0)
4742           sh_arg1=pokeShape(arg1)
4743           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4744           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4745           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4746           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4747           if not sh01==sh10:
4748               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)
4749           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4750    
4751       def getMyCode(self,argstrs,format="escript"):
4752          """
4753          returns a program code that can be used to evaluate the symbol.
4754    
4755          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4756          @type argstrs: C{list} of length 2 of C{str}.
4757          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4758          @type format: C{str}
4759          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4760          @rtype: C{str}
4761          @raise NotImplementedError: if the requested format is not available
4762          """
4763          if format=="escript" or format=="str" or format=="text":
4764             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4765          else:
4766             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4767    
4768       def substitute(self,argvals):
4769          """
4770          assigns new values to symbols in the definition of the symbol.
4771          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4772    
4773          @param argvals: new values assigned to symbols
4774          @type argvals: C{dict} with keywords of type L{Symbol}.
4775          @return: result of the substitution process. Operations are executed as much as possible.
4776          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4777          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4778          """
4779          if argvals.has_key(self):
4780             arg=argvals[self]
4781             if self.isAppropriateValue(arg):
4782                return arg
4783             else:
4784                raise TypeError,"%s: new value is not appropriate."%str(self)
4785          else:
4786             args=self.getSubstitutedArguments(argvals)
4787             return generalTensorTransposedProduct(args[0],args[1],args[2])
4788    
4789    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4790        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4791        # calculate the return shape:
4792        shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4793        shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4794        shape10=arg1.getShape()[arg1.getRank()-axis_offset:]
4795        shape1=arg1.getShape()[:arg1.getRank()-axis_offset]
4796        if not shape01==shape10:
4797            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)
4798    
4799        # whatr function space should be used? (this here is not good!)
4800        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4801        # create return value:
4802        out=escript.Data(0.,tuple(shape0+shape1),fs)
4803        #
4804        s0=[[]]
4805        for k in shape0:
4806              s=[]
4807              for j in s0:
4808                    for i in range(k): s.append(j+[slice(i,i)])
4809              s0=s
4810        s1=[[]]
4811        for k in shape1:
4812              s=[]
4813              for j in s1:
4814                    for i in range(k): s.append(j+[slice(i,i)])
4815              s1=s
4816        s01=[[]]
4817        for k in shape01:
4818              s=[]
4819              for j in s01:
4820                    for i in range(k): s.append(j+[slice(i,i)])
4821              s01=s
4822    
4823        for i0 in s0:
4824           for i1 in s1:
4825             s=escript.Scalar(0.,fs)
4826             for i01 in s01:
4827                s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i1+i01))
4828             out.__setitem__(tuple(i0+i1),s)
4829        return out
4830    
4831    
4832  #=========================================================  #=========================================================
4833  #   some little helpers  #  functions dealing with spatial dependency
4834  #=========================================================  #=========================================================
4835  def grad(arg,where=None):  def grad(arg,where=None):
4836      """      """
4837      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4838    
4839      @param arg:   Data object representing the function which gradient      If C{g} is the returned object, then
4840                    to be calculated.  
4841          - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4842          - if C{arg} is rank 1 C{g[i,s]} is the derivative of C{arg[i]} with respect to the C{s}-th spatial dimension.
4843          - if C{arg} is rank 2 C{g[i,j,s]} is the derivative of C{arg[i,j]} with respect to the C{s}-th spatial dimension.
4844          - if C{arg} is rank 3 C{g[i,j,k,s]} is the derivative of C{arg[i,j,k]} with respect to the C{s}-th spatial dimension.
4845    
4846        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4847        @type arg: L{escript.Data} or L{Symbol}
4848      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4849                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4850        @type where: C{None} or L{escript.FunctionSpace}
4851        @return: gradient of arg.
4852        @rtype: L{escript.Data} or L{Symbol}
4853      """      """
4854      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4855         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3617  def grad(arg,where=None): Line 4859  def grad(arg,where=None):
4859         else:         else:
4860            return arg._grad(where)            return arg._grad(where)
4861      else:      else:
4862        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4863    
4864    class Grad_Symbol(DependendSymbol):
4865       """
4866       L{Symbol} representing the result of the gradient operator
4867       """
4868       def __init__(self,arg,where=None):
4869          """
4870          initialization of gradient L{Symbol} with argument arg
4871          @param arg: argument of function
4872          @type arg: L{Symbol}.
4873          @param where: FunctionSpace in which the gradient will be calculated.
4874                      If not present or C{None} an appropriate default is used.
4875          @type where: C{None} or L{escript.FunctionSpace}
4876          """
4877          d=arg.getDim()
4878          if d==None:
4879             raise ValueError,"argument must have a spatial dimension"
4880          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4881    
4882       def getMyCode(self,argstrs,format="escript"):
4883          """
4884          returns a program code that can be used to evaluate the symbol.
4885    
4886          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4887          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4888          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4889          @type format: C{str}
4890          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4891          @rtype: C{str}
4892          @raise NotImplementedError: if the requested format is not available
4893          """
4894          if format=="escript" or format=="str"  or format=="text":
4895             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4896          else:
4897             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4898    
4899       def substitute(self,argvals):
4900          """
4901          assigns new values to symbols in the definition of the symbol.
4902          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4903    
4904          @param argvals: new values assigned to symbols
4905          @type argvals: C{dict} with keywords of type L{Symbol}.
4906          @return: result of the substitution process. Operations are executed as much as possible.
4907          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4908          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4909          """
4910          if argvals.has_key(self):
4911             arg=argvals[self]
4912             if self.isAppropriateValue(arg):
4913                return arg
4914             else:
4915                raise TypeError,"%s: new value is not appropriate."%str(self)
4916          else:
4917             arg=self.getSubstitutedArguments(argvals)
4918             return grad(arg[0],where=arg[1])
4919    
4920       def diff(self,arg):
4921          """
4922          differential of this object
4923    
4924          @param arg: the derivative is calculated with respect to arg
4925          @type arg: L{escript.Symbol}
4926          @return: derivative with respect to C{arg}
4927          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4928          """
4929          if arg==self:
4930             return identity(self.getShape())
4931          else:
4932             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4933    
4934  def integrate(arg,where=None):  def integrate(arg,where=None):
4935      """      """
4936      Return the integral if the function represented by Data object arg over      Return the integral of the function C{arg} over its domain. If C{where} is present C{arg} is interpolated to C{where}
4937      its domain.      before integration.
4938    
4939      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4940        @type arg: L{escript.Data} or L{Symbol}
4941      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4942                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4943        @type where: C{None} or L{escript.FunctionSpace}
4944        @return: integral of arg.
4945        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4946      """      """
4947      if isinstance(arg,numarray.NumArray):      if isinstance(arg,Symbol):
         if checkForZero(arg):  
            return arg  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,float):  
         if checkForZero(arg):  
            return arg  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,int):  
         if checkForZero(arg):  
            return float(arg)  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,Symbol):  
4948         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4949      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4950         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
# Line 3654  def integrate(arg,where=None): Line 4955  def integrate(arg,where=None):
4955      else:      else:
4956        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4957    
4958    class Integrate_Symbol(DependendSymbol):
4959       """
4960       L{Symbol} representing the result of the spatial integration operator
4961       """
4962       def __init__(self,arg,where=None):
4963          """
4964          initialization of integration L{Symbol} with argument arg
4965          @param arg: argument of the integration
4966          @type arg: L{Symbol}.
4967          @param where: FunctionSpace in which the integration will be calculated.
4968                      If not present or C{None} an appropriate default is used.
4969          @type where: C{None} or L{escript.FunctionSpace}
4970          """
4971          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4972    
4973       def getMyCode(self,argstrs,format="escript"):
4974          """
4975          returns a program code that can be used to evaluate the symbol.
4976    
4977          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4978          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4979          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4980          @type format: C{str}
4981          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4982          @rtype: C{str}
4983          @raise NotImplementedError: if the requested format is not available
4984          """
4985          if format=="escript" or format=="str"  or format=="text":
4986             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4987          else:
4988             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4989    
4990       def substitute(self,argvals):
4991          """
4992          assigns new values to symbols in the definition of the symbol.
4993          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4994    
4995          @param argvals: new values assigned to symbols
4996          @type argvals: C{dict} with keywords of type L{Symbol}.
4997          @return: result of the substitution process. Operations are executed as much as possible.
4998          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4999          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
5000          """
5001          if argvals.has_key(self):
5002             arg=argvals[self]
5003             if self.isAppropriateValue(arg):
5004                return arg
5005             else:
5006                raise TypeError,"%s: new value is not appropriate."%str(self)
5007          else:
5008             arg=self.getSubstitutedArguments(argvals)
5009             return integrate(arg[0],where=arg[1])
5010    
5011       def diff(self,arg):
5012          """
5013          differential of this object
5014    
5015          @param arg: the derivative is calculated with respect to arg
5016          @type arg: L{escript.Symbol}
5017          @return: derivative with respect to C{arg}
5018          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
5019          """
5020          if arg==self:
5021             return identity(self.getShape())
5022          else:
5023             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5024    
5025    
5026  def interpolate(arg,where):  def interpolate(arg,where):
5027      """      """
5028      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
5029    
5030      @param arg:    interpolant      @param arg: interpolant
5031      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
5032        @param where: FunctionSpace to be interpolated to
5033        @type where: L{escript.FunctionSpace}
5034        @return: interpolated argument
5035        @rtype: C{escript.Data} or L{Symbol}
5036      """      """
5037      if testForZero(arg):      if isinstance(arg,Symbol):
5038        return 0         return Interpolate_Symbol(arg,where)
     elif isinstance(arg,Symbol):  
        return Interpolated_Symbol(arg,where)  
5039      else:      else:
5040         return escript.Data(arg,where)         return escript.Data(arg,where)
5041    
5042  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
5043      """     """
5044      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
5045       """
5046       def __init__(self,arg,where):
5047          """
5048          initialization of interpolation L{Symbol} with argument arg
5049          @param arg: argument of the interpolation
5050          @type arg: L{Symbol}.
5051          @param where: FunctionSpace into which the argument is interpolated.
5052          @type where: L{escript.FunctionSpace}
5053          """
5054          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
5055    
5056      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
5057                    be calculated.        """
5058      @param where: FunctionSpace in which the gradient will be calculated.        returns a program code that can be used to evaluate the symbol.
                   If not present or C{None} an appropriate default is used.  
     """  
     g=grad(arg,where)  
     return trace(g,axis0=g.getRank()-2,axis1=g.getRank()-1)  
5059    
5060  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
5061      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
5062      Returns the jump of arg across a continuity.        @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
5063          @type format: C{str}
5064          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
5065          @rtype: C{str}
5066          @raise NotImplementedError: if the requested format is not available
5067          """
5068          if format=="escript" or format=="str"  or format=="text":
5069             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
5070          else:
5071             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
5072    
5073      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
5074                    to be calculated.        """
5075      """        assigns new values to symbols in the definition of the symbol.
5076      d=arg.getDomain()        The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
     return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())  
5077    
5078  #=============================        @param argvals: new values assigned to symbols
5079  #        @type argvals: C{dict} with keywords of type L{Symbol}.
5080  # wrapper for various functions: if the argument has attribute the function name        @return: result of the substitution process. Operations are executed as much as possible.
5081  # as an argument it calls the corresponding methods. Otherwise the corresponding        @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
5082  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
5083          """
5084          if argvals.has_key(self):
5085             arg=argvals[self]
5086             if self.isAppropriateValue(arg):
5087                return arg
5088             else:
5089                raise TypeError,"%s: new value is not appropriate."%str(self)
5090          else:
5091             arg=self.getSubstitutedArguments(argvals)
5092             return interpolate(arg[0],where=arg[1])
5093    
5094       def diff(self,arg):
5095          """
5096          differential of this object
5097    
5098          @param arg: the derivative is calculated with respect to arg
5099          @type arg: L{escript.Symbol}
5100          @return: derivative with respect to C{arg}
5101          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
5102          """
5103          if arg==self:
5104             return identity(self.getShape())
5105          else:
5106             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5107    
 # functions involving the underlying Domain:  
5108    
5109  def transpose(arg,axis=None):  def div(arg,where=None):
5110      """      """
5111      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
5112    
5113      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
5114        @type arg: L{escript.Data} or L{Symbol}
5115        @param where: FunctionSpace in which the divergence will be calculated.
5116                      If not present or C{None} an appropriate default is used.
5117        @type where: C{None} or L{escript.FunctionSpace}
5118        @return: divergence of arg.
5119        @rtype: L{escript.Data} or L{Symbol}
5120      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
5121      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5122         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
5123      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
5124         # hack for transpose          dim=arg.getDomain().getDim()
        r=arg.getRank()  
        if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"  
        s=arg.getShape()  
        out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())  
        for i in range(s[0]):  
           for j in range(s[1]):  
              out[j,i]=arg[i,j]  
        return out  
        # end hack for transpose  
        return arg.transpose(axis)  
5125      else:      else:
5126         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
5127        if not arg.getShape()==(dim,):
5128          raise ValueError,"div: expected shape is (%s,)"%dim
5129        return trace(grad(arg,where))
5130    
5131  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
5132      """      """
5133      Return      returns the jump of arg across the continuity of the domain
5134    
5135      @param arg:      @param arg: argument
5136        @type arg: L{escript.Data} or L{Symbol}
5137        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5138                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5139        @type domain: C{None} or L{escript.Domain}
5140        @return: jump of arg
5141        @rtype: L{escript.Data} or L{Symbol}
5142      """      """
5143      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
5144         s=list(arg.getShape())              return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
        s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])  
        return Trace_Symbol(arg,axis0=axis0,axis1=axis1)  
     elif isinstance(arg,escript.Data):  
        # hack for trace  
        s=arg.getShape()  
        if s[axis0]!=s[axis1]:  
            raise ValueError,"illegal axis in trace"  
        out=escript.Scalar(0.,arg.getFunctionSpace())  
        for i in range(s[axis0]):  
           out+=arg[i,i]  
        return out  
        # end hack for trace  
     else:  
        return numarray.trace(arg,axis0=axis0,axis1=axis1)  
5145    
5146    def L2(arg):
5147        """
5148        returns the L2 norm of arg at where
5149        
5150        @param arg: function which L2 to be calculated.
5151        @type arg: L{escript.Data} or L{Symbol}
5152        @return: L2 norm of arg.
5153        @rtype: L{float} or L{Symbol}
5154        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5155        """
5156        return sqrt(integrate(inner(arg,arg)))
5157    #=============================
5158    #
5159    
5160  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5161      """      """
5162      resorts the component of arg according to index      resorts the component of arg according to index
5163    
5164      """      """
5165      pass      raise NotImplementedError
5166  #  #
5167  # $Log: util.py,v $  # $Log: util.py,v $
5168  # Revision 1.14.2.16  2005/10/19 06:09:57  gross  # Revision 1.14.2.16  2005/10/19 06:09:57  gross

Legend:
Removed from v.341  
changed lines
  Added in v.804

  ViewVC Help
Powered by ViewVC 1.1.26