/[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 396 by gross, Wed Dec 21 05:08:25 2005 UTC revision 775 by ksteube, Mon Jul 10 04:00:08 2006 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2    
3  """  """
4  Utility functions for escript  Utility functions for escript
5    
 @remark:  This module is under construction and is still tested!!!  
   
6  @var __author__: name of author  @var __author__: name of author
7  @var __licence__: licence agreement  @var __copyright__: copyrights
8    @var __license__: licence agreement
9  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
10  @var __version__: version  @var __version__: version
11  @var __date__: date of the version  @var __date__: date of the version
12  """  """
13                                                                                                                                                                                                                                                                                                                                                                                                            
14  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
15  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                        http://www.access.edu.au
17                    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19                 http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
21  __version__="$Revision$"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
# Line 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 clip(arg,minval,maxval)  
   
 # 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 69  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 97  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 126  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 144  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 164  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 179  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 192  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 250  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 339  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 361  def testForZero(arg): Line 330  def testForZero(arg):
330    
331      @param arg: a given object      @param arg: a given object
332      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}
333      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
334      @rtype : C{bool}      @rtype: C{bool}
335      """      """
336      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
337         return not Lsup(arg)>0.         return not Lsup(arg)>0.
# Line 395  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 421  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 465  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
438    
439      If shape is not given the shape "largest" shape of args is used.      @param arg0: a given object
440        @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
441      @param args: a given ob      @param arg1: a given object
442      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
443      @return: True if the argument is identical to zero.      @return: arg0 and arg1 where copies are returned when the shape has to be changed.
444      @rtype: C{list} of C{int}      @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 500  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 544  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 553  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 577  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 606  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 696  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 704  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 713  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 730  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 741  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 752  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 763  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 774  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 785  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 796  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 807  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 818  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 829  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 884  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 912  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        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
988        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
989        return out        return out
990     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
991        return arg._wherePositive()        return arg._wherePositive()
# Line 958  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 994  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        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1070        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1071        return out        return out
1072     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1073        return arg._whereNegative()        return arg._whereNegative()
# Line 1040  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 1076  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        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1152        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1153        return out        return out
1154     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1155        return arg._whereNonNegative()        return arg._whereNonNegative()
# Line 1106  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        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1182        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1183        return out        return out
1184     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1185        return arg._whereNonPositive()        return arg._whereNonPositive()
# Line 1138  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        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1214        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1215        return out        return out
1216     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1217        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1218     elif isinstance(arg,float):     elif isinstance(arg,float):
1219        if abs(arg)<=tol:        if abs(arg)<=tol:
1220          return 1.          return 1.
# Line 1187  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 1221  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        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1294        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1295        return out        return out
1296     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1297        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1298     elif isinstance(arg,float):     elif isinstance(arg,float):
1299        if abs(arg)>tol:        if abs(arg)>tol:
1300          return 1.          return 1.
# Line 1254  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 1292  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 1344  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 1382  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 1434  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 1472  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 1524  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 1562  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 1614  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 1652  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 1704  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 1742  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 1794  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 1832  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 1884  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 1922  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 1974  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 2012  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 2064  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 2102  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 2154  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 2192  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 2244  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 2282  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 2334  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 2372  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 2424  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 2462  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 2514  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 2552  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 2604  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 2652  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 2704  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 2745  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 2781  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 2822  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 2858  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: 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                      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,"trace: rank of argument must be greater than 1"
2943          if axis_offset<0 or axis_offset>len(sh)-2:
2944            raise ValueError,"trace: axis_offset must be between 0 and %s"%len(sh)-2
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,"trace: 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          return escript_trace(arg,axis_offset)
2960       elif isinstance(arg,float):
2961          raise TypeError,"trace: illegal argument type float."
2962       elif isinstance(arg,int):
2963          raise TypeError,"trace: illegal argument type int."
2964       elif isinstance(arg,Symbol):
2965          return Trace_Symbol(arg,axis_offset)
2966       else:
2967          raise TypeError,"trace: Unknown argument type."
2968    
2969    def escript_trace(arg,axis_offset):
2970          "arg is a Data object"
2971          if arg.getRank()<2:
2972            raise ValueError,"escript_trace: rank of argument must be greater than 1"
2973          if axis_offset<0 or axis_offset>arg.getRank()-2:
2974            raise ValueError,"escript_trace: axis_offset must be between 0 and %s"%arg.getRank()-2
2975          s=list(arg.getShape())        
2976          if not s[axis_offset] == s[axis_offset+1]:
2977            raise ValueError,"escript_trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2978          out=arg._matrixtrace(axis_offset)
2979          return out
2980    class Trace_Symbol(DependendSymbol):
2981       """
2982       L{Symbol} representing the result of the trace function
2983       """
2984       def __init__(self,arg,axis_offset=0):
2985          """
2986          initialization of trace L{Symbol} with argument arg
2987          @param arg: argument of function
2988          @type arg: L{Symbol}.
2989          @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
2990                      axis_offset and axis_offset+1 must be equal.
2991          @type axis_offset: C{int}
2992          """
2993          if arg.getRank()<2:
2994            raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"
2995          if axis_offset<0 or axis_offset>arg.getRank()-2:
2996            raise ValueError,"Trace_Symbol: axis_offset must be between 0 and %s"%arg.getRank()-2
2997          s=list(arg.getShape())        
2998          if not s[axis_offset] == s[axis_offset+1]:
2999            raise ValueError,"Trace_Symbol: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3000          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3001    
3002       def getMyCode(self,argstrs,format="escript"):
3003          """
3004          returns a program code that can be used to evaluate the symbol.
3005    
3006          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3007          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3008          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3009          @type format: C{str}
3010          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3011          @rtype: C{str}
3012          @raise NotImplementedError: if the requested format is not available
3013          """
3014          if format=="escript" or format=="str"  or format=="text":
3015             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3016          else:
3017             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3018    
3019       def substitute(self,argvals):
3020          """
3021          assigns new values to symbols in the definition of the symbol.
3022          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3023    
3024          @param argvals: new values assigned to symbols
3025          @type argvals: C{dict} with keywords of type L{Symbol}.
3026          @return: result of the substitution process. Operations are executed as much as possible.
3027          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3028          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3029          """
3030          if argvals.has_key(self):
3031             arg=argvals[self]
3032             if self.isAppropriateValue(arg):
3033                return arg
3034             else:
3035                raise TypeError,"%s: new value is not appropriate."%str(self)
3036          else:
3037             arg=self.getSubstitutedArguments(argvals)
3038             return trace(arg[0],axis_offset=arg[1])
3039    
3040       def diff(self,arg):
3041          """
3042          differential of this object
3043    
3044          @param arg: the derivative is calculated with respect to arg
3045          @type arg: L{escript.Symbol}
3046          @return: derivative with respect to C{arg}
3047          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3048          """
3049          if arg==self:
3050             return identity(self.getShape())
3051          else:
3052             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3053    
3054    def transpose(arg,axis_offset=None):
3055       """
3056       returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.
3057    
3058       @param arg: argument
3059       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3060       @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3061                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3062       @type axis_offset: C{int}
3063       @return: transpose of arg
3064       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3065       """
3066       if isinstance(arg,numarray.NumArray):
3067          if axis_offset==None: axis_offset=int(arg.rank/2)
3068          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3069       elif isinstance(arg,escript.Data):
3070          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3071          return escript_transpose(arg,axis_offset)
3072       elif isinstance(arg,float):
3073          if not ( axis_offset==0 or axis_offset==None):
3074            raise ValueError,"transpose: axis_offset must be 0 for float argument"
3075          return arg
3076       elif isinstance(arg,int):
3077          if not ( axis_offset==0 or axis_offset==None):
3078            raise ValueError,"transpose: axis_offset must be 0 for int argument"
3079          return float(arg)
3080       elif isinstance(arg,Symbol):
3081          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3082          return Transpose_Symbol(arg,axis_offset)
3083       else:
3084          raise TypeError,"transpose: Unknown argument type."
3085    
3086    def escript_transpose(arg,axis_offset):
3087          "arg is a Data object"
3088          r=arg.getRank()
3089          if axis_offset<0 or axis_offset>r:
3090            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3091          out=arg._transpose(axis_offset)
3092          return out
3093    class Transpose_Symbol(DependendSymbol):
3094       """
3095       L{Symbol} representing the result of the transpose function
3096       """
3097       def __init__(self,arg,axis_offset=None):
3098          """
3099          initialization of transpose L{Symbol} with argument arg
3100    
3101          @param arg: argument of function
3102          @type arg: L{Symbol}.
3103          @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3104                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3105          @type axis_offset: C{int}
3106          """
3107          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3108          if axis_offset<0 or axis_offset>arg.getRank():
3109            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3110          s=arg.getShape()
3111          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3112    
3113       def getMyCode(self,argstrs,format="escript"):
3114          """
3115          returns a program code that can be used to evaluate the symbol.
3116    
3117          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3118          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3119          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3120          @type format: C{str}
3121          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3122          @rtype: C{str}
3123          @raise NotImplementedError: if the requested format is not available
3124          """
3125          if format=="escript" or format=="str"  or format=="text":
3126             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3127          else:
3128             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3129    
3130       def substitute(self,argvals):
3131          """
3132          assigns new values to symbols in the definition of the symbol.
3133          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3134    
3135          @param argvals: new values assigned to symbols
3136          @type argvals: C{dict} with keywords of type L{Symbol}.
3137          @return: result of the substitution process. Operations are executed as much as possible.
3138          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3139          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3140          """
3141          if argvals.has_key(self):
3142             arg=argvals[self]
3143             if self.isAppropriateValue(arg):
3144                return arg
3145             else:
3146                raise TypeError,"%s: new value is not appropriate."%str(self)
3147          else:
3148             arg=self.getSubstitutedArguments(argvals)
3149             return transpose(arg[0],axis_offset=arg[1])
3150    
3151       def diff(self,arg):
3152          """
3153          differential of this object
3154    
3155          @param arg: the derivative is calculated with respect to arg
3156          @type arg: L{escript.Symbol}
3157          @return: derivative with respect to C{arg}
3158          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3159          """
3160          if arg==self:
3161             return identity(self.getShape())
3162          else:
3163             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3164    def symmetric(arg):
3165        """
3166        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3167    
3168        @param arg: square matrix. Must have rank 2 or 4 and be square.
3169        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3170        @return: symmetric part of arg
3171        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3172        """
3173        if isinstance(arg,numarray.NumArray):
3174          if arg.rank==2:
3175            if not (arg.shape[0]==arg.shape[1]):
3176               raise ValueError,"symmetric: argument must be square."
3177          elif arg.rank==4:
3178            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3179               raise ValueError,"symmetric: argument must be square."
3180          else:
3181            raise ValueError,"symmetric: rank 2 or 4 is required."
3182          return (arg+transpose(arg))/2
3183        elif isinstance(arg,escript.Data):
3184          return escript_symmetric(arg)
3185        elif isinstance(arg,float):
3186          return arg
3187        elif isinstance(arg,int):
3188          return float(arg)
3189        elif isinstance(arg,Symbol):
3190          if arg.getRank()==2:
3191            if not (arg.getShape()[0]==arg.getShape()[1]):
3192               raise ValueError,"symmetric: argument must be square."
3193          elif arg.getRank()==4:
3194            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3195               raise ValueError,"symmetric: argument must be square."
3196          else:
3197            raise ValueError,"symmetric: rank 2 or 4 is required."
3198          return (arg+transpose(arg))/2
3199        else:
3200          raise TypeError,"symmetric: Unknown argument type."
3201    
3202    def escript_symmetric(arg):
3203          if arg.getRank()==2:
3204            if not (arg.getShape()[0]==arg.getShape()[1]):
3205               raise ValueError,"escript_symmetric: argument must be square."
3206            out=arg._symmetric()
3207          elif arg.getRank()==4:
3208            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3209               raise ValueError,"escript_symmetric: argument must be square."
3210            out=arg._symmetric()
3211          else:
3212            raise ValueError,"escript_symmetric: rank 2 or 4 is required."
3213          return out
3214    
3215    def nonsymmetric(arg):
3216        """
3217        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3218    
3219        @param arg: square matrix. Must have rank 2 or 4 and be square.
3220        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3221        @return: nonsymmetric part of arg
3222        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3223        """
3224        if isinstance(arg,numarray.NumArray):
3225          if arg.rank==2:
3226            if not (arg.shape[0]==arg.shape[1]):
3227               raise ValueError,"nonsymmetric: argument must be square."
3228          elif arg.rank==4:
3229            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3230               raise ValueError,"nonsymmetric: argument must be square."
3231          else:
3232            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3233          return (arg-transpose(arg))/2
3234        elif isinstance(arg,escript.Data):
3235          return escript_nonsymmetric(arg)
3236        elif isinstance(arg,float):
3237          return arg
3238        elif isinstance(arg,int):
3239          return float(arg)
3240        elif isinstance(arg,Symbol):
3241          if arg.getRank()==2:
3242            if not (arg.getShape()[0]==arg.getShape()[1]):
3243               raise ValueError,"nonsymmetric: argument must be square."
3244          elif arg.getRank()==4:
3245            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3246               raise ValueError,"nonsymmetric: argument must be square."
3247          else:
3248            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3249          return (arg-transpose(arg))/2
3250        else:
3251          raise TypeError,"nonsymmetric: Unknown argument type."
3252    
3253    def escript_nonsymmetric(arg): # this should be implemented in c++
3254          if arg.getRank()==2:
3255            if not (arg.getShape()[0]==arg.getShape()[1]):
3256               raise ValueError,"escript_nonsymmetric: argument must be square."
3257            out=arg._nonsymmetric()
3258          elif arg.getRank()==4:
3259            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3260               raise ValueError,"escript_nonsymmetric: argument must be square."
3261            out=arg._nonsymmetric()
3262          else:
3263            raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."
3264          return out
3265    
3266    
3267    def inverse(arg):
3268        """
3269        returns the inverse of the square matrix arg.
3270    
3271        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3272        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3273        @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3274        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3275        @note: for L{escript.Data} objects the dimension is restricted to 3.
3276        """
3277        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3278        if isinstance(arg,numarray.NumArray):
3279          return numarray.linear_algebra.inverse(arg)
3280        elif isinstance(arg,escript.Data):
3281          return escript_inverse(arg)
3282        elif isinstance(arg,float):
3283          return 1./arg
3284        elif isinstance(arg,int):
3285          return 1./float(arg)
3286        elif isinstance(arg,Symbol):
3287          return Inverse_Symbol(arg)
3288        else:
3289          raise TypeError,"inverse: Unknown argument type."
3290    
3291    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3292          "arg is a Data objects!!!"
3293          if not arg.getRank()==2:
3294            raise ValueError,"escript_inverse: argument must have rank 2"
3295          s=arg.getShape()      
3296          if not s[0] == s[1]:
3297            raise ValueError,"escript_inverse: argument must be a square matrix."
3298          out=escript.Data(0.,s,arg.getFunctionSpace())
3299          if s[0]==1:
3300              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3301                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3302              out[0,0]=1./arg[0,0]
3303          elif s[0]==2:
3304              A11=arg[0,0]
3305              A12=arg[0,1]
3306              A21=arg[1,0]
3307              A22=arg[1,1]
3308              D = A11*A22-A12*A21
3309              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3310                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3311              D=1./D
3312              out[0,0]= A22*D
3313              out[1,0]=-A21*D
3314              out[0,1]=-A12*D
3315              out[1,1]= A11*D
3316          elif s[0]==3:
3317              A11=arg[0,0]
3318              A21=arg[1,0]
3319              A31=arg[2,0]
3320              A12=arg[0,1]
3321              A22=arg[1,1]
3322              A32=arg[2,1]
3323              A13=arg[0,2]
3324              A23=arg[1,2]
3325              A33=arg[2,2]
3326              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3327              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3328                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3329              D=1./D
3330              out[0,0]=(A22*A33-A23*A32)*D
3331              out[1,0]=(A31*A23-A21*A33)*D
3332              out[2,0]=(A21*A32-A31*A22)*D
3333              out[0,1]=(A13*A32-A12*A33)*D
3334              out[1,1]=(A11*A33-A31*A13)*D
3335              out[2,1]=(A12*A31-A11*A32)*D
3336              out[0,2]=(A12*A23-A13*A22)*D
3337              out[1,2]=(A13*A21-A11*A23)*D
3338              out[2,2]=(A11*A22-A12*A21)*D
3339          else:
3340             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3341          return out
3342    
3343    class Inverse_Symbol(DependendSymbol):
3344       """
3345       L{Symbol} representing the result of the inverse function
3346       """
3347       def __init__(self,arg):
3348          """
3349          initialization of inverse L{Symbol} with argument arg
3350          @param arg: argument of function
3351          @type arg: L{Symbol}.
3352          """
3353          if not arg.getRank()==2:
3354            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3355          s=arg.getShape()
3356          if not s[0] == s[1]:
3357            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3358          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3359    
3360       def getMyCode(self,argstrs,format="escript"):
3361          """
3362          returns a program code that can be used to evaluate the symbol.
3363    
3364          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3365          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3366          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3367          @type format: C{str}
3368          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3369          @rtype: C{str}
3370          @raise NotImplementedError: if the requested format is not available
3371          """
3372          if format=="escript" or format=="str"  or format=="text":
3373             return "inverse(%s)"%argstrs[0]
3374          else:
3375             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3376    
3377       def substitute(self,argvals):
3378          """
3379          assigns new values to symbols in the definition of the symbol.
3380          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3381    
3382          @param argvals: new values assigned to symbols
3383          @type argvals: C{dict} with keywords of type L{Symbol}.
3384          @return: result of the substitution process. Operations are executed as much as possible.
3385          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3386          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3387          """
3388          if argvals.has_key(self):
3389             arg=argvals[self]
3390             if self.isAppropriateValue(arg):
3391                return arg
3392             else:
3393                raise TypeError,"%s: new value is not appropriate."%str(self)
3394          else:
3395             arg=self.getSubstitutedArguments(argvals)
3396             return inverse(arg[0])
3397    
3398       def diff(self,arg):
3399          """
3400          differential of this object
3401    
3402          @param arg: the derivative is calculated with respect to arg
3403          @type arg: L{escript.Symbol}
3404          @return: derivative with respect to C{arg}
3405          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3406          """
3407          if arg==self:
3408             return identity(self.getShape())
3409          else:
3410             return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)
3411    
3412    def eigenvalues(arg):
3413        """
3414        returns the eigenvalues of the square matrix arg.
3415    
3416        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3417                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3418        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3419        @return: the eigenvalues in increasing order.
3420        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3421        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3422        """
3423        if isinstance(arg,numarray.NumArray):
3424          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3425          out.sort()
3426          return out
3427        elif isinstance(arg,escript.Data):
3428          return arg._eigenvalues()
3429        elif isinstance(arg,Symbol):
3430          if not arg.getRank()==2:
3431            raise ValueError,"eigenvalues: argument must have rank 2"
3432          s=arg.getShape()      
3433          if not s[0] == s[1]:
3434            raise ValueError,"eigenvalues: argument must be a square matrix."
3435          if s[0]==1:
3436              return arg[0]
3437          elif s[0]==2:
3438              arg1=symmetric(arg)
3439              A11=arg1[0,0]
3440              A12=arg1[0,1]
3441              A22=arg1[1,1]
3442              trA=(A11+A22)/2.
3443              A11-=trA
3444              A22-=trA
3445              s=sqrt(A12**2-A11*A22)
3446              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3447          elif s[0]==3:
3448              arg1=symmetric(arg)
3449              A11=arg1[0,0]
3450              A12=arg1[0,1]
3451              A22=arg1[1,1]
3452              A13=arg1[0,2]
3453              A23=arg1[1,2]
3454              A33=arg1[2,2]
3455              trA=(A11+A22+A33)/3.
3456              A11-=trA
3457              A22-=trA
3458              A33-=trA
3459              A13_2=A13**2
3460              A23_2=A23**2
3461              A12_2=A12**2
3462              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3463              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3464              sq_p=sqrt(p/3.)
3465              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
3466              sq_p*=2.
3467              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3468               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3469               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3470              return trA+sq_p*f
3471          else:
3472             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3473        elif isinstance(arg,float):
3474          return arg
3475        elif isinstance(arg,int):
3476          return float(arg)
3477        else:
3478          raise TypeError,"eigenvalues: Unknown argument type."
3479    
3480    def eigenvalues_and_eigenvectors(arg):
3481        """
3482        returns the eigenvalues and eigenvectors of the square matrix arg.
3483    
3484        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3485                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3486        @type arg: L{escript.Data}
3487        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3488                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3489                 the eigenvector coresponding to the i-th eigenvalue.
3490        @rtype: L{tuple} of L{escript.Data}.
3491        @note: The dimension is restricted to 3.
3492        """
3493        if isinstance(arg,numarray.NumArray):
3494          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3495        elif isinstance(arg,escript.Data):
3496          return arg._eigenvalues_and_eigenvectors()
3497        elif isinstance(arg,Symbol):
3498          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3499        elif isinstance(arg,float):
3500          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3501        elif isinstance(arg,int):
3502          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3503        else:
3504          raise TypeError,"eigenvalues: Unknown argument type."
3505  #=======================================================  #=======================================================
3506  #  Binary operations:  #  Binary operations:
3507  #=======================================================  #=======================================================
# Line 2921  class Add_Symbol(DependendSymbol): Line 3561  class Add_Symbol(DependendSymbol):
3561        @type format: C{str}        @type format: C{str}
3562        @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.
3563        @rtype: C{str}        @rtype: C{str}
3564        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3565        """        """
3566        if format=="str" or format=="text":        if format=="str" or format=="text":
3567           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 2980  def mult(arg0,arg1): Line 3620  def mult(arg0,arg1):
3620         """         """
3621         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3622         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3623            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3624         else:         else:
3625            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3626                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3020  class Mult_Symbol(DependendSymbol): Line 3660  class Mult_Symbol(DependendSymbol):
3660        @type format: C{str}        @type format: C{str}
3661        @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.
3662        @rtype: C{str}        @rtype: C{str}
3663        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3664        """        """
3665        if format=="str" or format=="text":        if format=="str" or format=="text":
3666           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3080  def quotient(arg0,arg1): Line 3720  def quotient(arg0,arg1):
3720         """         """
3721         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3722         if testForZero(args[0]):         if testForZero(args[0]):
3723            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3724         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3725            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3726               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3125  class Quotient_Symbol(DependendSymbol): Line 3765  class Quotient_Symbol(DependendSymbol):
3765        @type format: C{str}        @type format: C{str}
3766        @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.
3767        @rtype: C{str}        @rtype: C{str}
3768        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3769        """        """
3770        if format=="str" or format=="text":        if format=="str" or format=="text":
3771           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3186  def power(arg0,arg1): Line 3826  def power(arg0,arg1):
3826         """         """
3827         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3828         if testForZero(args[0]):         if testForZero(args[0]):
3829            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3830         elif testForZero(args[1]):         elif testForZero(args[1]):
3831            return numarray.ones(args[0],numarray.Float)            return numarray.ones(pokeShape(args[1]),numarray.Float64)
3832         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3833            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3834         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 3229  class Power_Symbol(DependendSymbol): Line 3869  class Power_Symbol(DependendSymbol):
3869        @type format: C{str}        @type format: C{str}
3870        @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.
3871        @rtype: C{str}        @rtype: C{str}
3872        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3873        """        """
3874        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
3875           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3318  def clip(arg,minval=0.,maxval=1.): Line 3958  def clip(arg,minval=0.,maxval=1.):
3958      @param arg: argument      @param arg: argument
3959      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
3960      @param minval: lower range      @param minval: lower range
3961      @type arg: C{float}      @type minval: C{float}
3962      @param maxval: uper range      @param maxval: upper range
3963      @type arg: C{float}      @type maxval: C{float}
3964      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
3965               less then maxval are unchanged.               less then maxval are unchanged.
3966      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
3967        @raise ValueError: if minval>maxval
3968      """      """
3969      if minval>maxval:      if minval>maxval:
3970         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
# Line 3344  def inner(arg0,arg1): Line 3985  def inner(arg0,arg1):
3985      @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}
3986      @param arg1: second argument      @param arg1: second argument
3987      @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}
3988      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
3989      @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
3990      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
3991      """      """
# Line 3352  def inner(arg0,arg1): Line 3993  def inner(arg0,arg1):
3993      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
3994      if not sh0==sh1:      if not sh0==sh1:
3995          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
3996      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
3997    
3998  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
3999      """      """
# Line 3360  def matrixmult(arg0,arg1): Line 4001  def matrixmult(arg0,arg1):
4001    
4002      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4003    
4004            or      or
4005    
4006      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4007    
# Line 3380  def matrixmult(arg0,arg1): Line 4021  def matrixmult(arg0,arg1):
4021          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4022      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4023          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4024      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4025    
4026  def outer(arg0,arg1):  def outer(arg0,arg1):
4027      """      """
# Line 3388  def outer(arg0,arg1): Line 4029  def outer(arg0,arg1):
4029    
4030      out[t,s]=arg0[t]*arg1[s]      out[t,s]=arg0[t]*arg1[s]
4031    
4032      where s runs through arg0.Shape      where
4033            t runs through arg1.Shape  
4034            - s runs through arg0.Shape
4035            - t runs through arg1.Shape
4036    
4037      @param arg0: first argument      @param arg0: first argument
4038      @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}
# Line 3398  def outer(arg0,arg1): Line 4041  def outer(arg0,arg1):
4041      @return: the outer product of arg0 and arg1 at each data point      @return: the outer product of arg0 and arg1 at each data point
4042      @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
4043      """      """
4044      return generalTensorProduct(arg0,arg1,offset=0)      return generalTensorProduct(arg0,arg1,axis_offset=0)
4045    
4046    
4047  def tensormult(arg0,arg1):  def tensormult(arg0,arg1):
4048      """      """
4049      the tensor product of the two argument:      the tensor product of the two argument:
   
4050            
4051      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4052    
4053      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4054    
4055                   or      or
4056    
4057      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4058    
# Line 3419  def tensormult(arg0,arg1): Line 4061  def tensormult(arg0,arg1):
4061    
4062      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]
4063                                
4064                   or      or
4065    
4066      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]
4067    
4068                   or      or
4069    
4070      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]
4071    
# Line 3440  def tensormult(arg0,arg1): Line 4082  def tensormult(arg0,arg1):
4082      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
4083      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4084      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4085         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4086      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):
4087         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4088      else:      else:
4089          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensormult: first argument must have rank 2 or 4"
4090    
4091  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4092      """      """
4093      generalized tensor product      generalized tensor product
4094    
4095      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4096    
4097      where s runs through arg0.Shape[:arg0.Rank-offset]      where
4098            r runs trough arg0.Shape[:offset]  
4099            t runs through arg1.Shape[offset:]          - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4100            - r runs trough arg0.Shape[:axis_offset]
4101            - t runs through arg1.Shape[axis_offset:]
4102    
4103      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 length of arg1 must match and  
4104      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 shape of arg1.
# Line 3471  def generalTensorProduct(arg0,arg1,offse Line 4115  def generalTensorProduct(arg0,arg1,offse
4115      # 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
4116      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4117         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4118             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4119         else:         else:
4120             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4121                 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,"generalTensorProduct: dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4122             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4123             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4124             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4125             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4126             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4127             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4128             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4129             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4130             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4131             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4132             for i0 in range(d0):             for i0 in range(d0):
4133                      for i1 in range(d1):                      for i1 in range(d1):
4134                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4135             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4136             return out             return out
4137      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4138         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4139             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4140         else:         else:
4141             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)
4142      else:            else:      
4143         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4144                                    
4145  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4146     """     """
4147     Symbol representing the quotient of two arguments.     Symbol representing the quotient of two arguments.
4148     """     """
4149     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4150         """         """
4151         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the quotient of two arguments
4152    
# Line 3515  class GeneralTensorProduct_Symbol(Depend Line 4159  class GeneralTensorProduct_Symbol(Depend
4159         """         """
4160         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
4161         sh_arg1=pokeShape(arg1)         sh_arg1=pokeShape(arg1)
4162         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4163         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4164         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4165         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4166         if not sh01==sh10:         if not sh01==sh10:
4167             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)
4168         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])
4169    
4170     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4171        """        """
# Line 3533  class GeneralTensorProduct_Symbol(Depend Line 4177  class GeneralTensorProduct_Symbol(Depend
4177        @type format: C{str}        @type format: C{str}
4178        @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.
4179        @rtype: C{str}        @rtype: C{str}
4180        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4181        """        """
4182        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4183           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])
4184        else:        else:
4185           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)
4186    
# Line 3561  class GeneralTensorProduct_Symbol(Depend Line 4205  class GeneralTensorProduct_Symbol(Depend
4205           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4206           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4207    
4208  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct
4209      "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!!!"
4210      # calculate the return shape:      # calculate the return shape:
4211      shape0=arg0.getShape()[:arg0.getRank()-offset]      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4212      shape01=arg0.getShape()[arg0.getRank()-offset:]      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4213      shape10=arg1.getShape()[:offset]      shape10=arg1.getShape()[:axis_offset]
4214      shape1=arg1.getShape()[offset:]      shape1=arg1.getShape()[axis_offset:]
4215      if not shape01==shape10:      if not shape01==shape10:
4216          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)
4217    
4218        # whatr function space should be used? (this here is not good!)
4219        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4220      # create return value:      # create return value:
4221      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())      out=escript.Data(0.,tuple(shape0+shape1),fs)
4222      #      #
4223      s0=[[]]      s0=[[]]
4224      for k in shape0:      for k in shape0:
# Line 3595  def escript_generalTensorProduct(arg0,ar Line 4241  def escript_generalTensorProduct(arg0,ar
4241    
4242      for i0 in s0:      for i0 in s0:
4243         for i1 in s1:         for i1 in s1:
4244           s=escript.Scalar(0.,arg0.getFunctionSpace())           s=escript.Scalar(0.,fs)
4245           for i01 in s01:           for i01 in s01:
4246              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))
4247           out.__setitem__(tuple(i0+i1),s)           out.__setitem__(tuple(i0+i1),s)
4248      return out      return out
4249    
4250    
4251  #=========================================================  #=========================================================
4252  #   some little helpers  #  functions dealing with spatial dependency
4253  #=========================================================  #=========================================================
4254  def grad(arg,where=None):  def grad(arg,where=None):
4255      """      """
4256      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4257    
4258      @param arg:   Data object representing the function which gradient      If C{g} is the returned object, then
4259                    to be calculated.  
4260          - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4261          - 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.
4262          - 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.
4263          - 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.
4264    
4265        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4266        @type arg: L{escript.Data} or L{Symbol}
4267      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4268                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4269        @type where: C{None} or L{escript.FunctionSpace}
4270        @return: gradient of arg.
4271        @rtype: L{escript.Data} or L{Symbol}
4272      """      """
4273      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4274         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3621  def grad(arg,where=None): Line 4278  def grad(arg,where=None):
4278         else:         else:
4279            return arg._grad(where)            return arg._grad(where)
4280      else:      else:
4281        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4282    
4283    class Grad_Symbol(DependendSymbol):
4284       """
4285       L{Symbol} representing the result of the gradient operator
4286       """
4287       def __init__(self,arg,where=None):
4288          """
4289          initialization of gradient L{Symbol} with argument arg
4290          @param arg: argument of function
4291          @type arg: L{Symbol}.
4292          @param where: FunctionSpace in which the gradient will be calculated.
4293                      If not present or C{None} an appropriate default is used.
4294          @type where: C{None} or L{escript.FunctionSpace}
4295          """
4296          d=arg.getDim()
4297          if d==None:
4298             raise ValueError,"argument must have a spatial dimension"
4299          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4300    
4301       def getMyCode(self,argstrs,format="escript"):
4302          """
4303          returns a program code that can be used to evaluate the symbol.
4304    
4305          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4306          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4307          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4308          @type format: C{str}
4309          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4310          @rtype: C{str}
4311          @raise NotImplementedError: if the requested format is not available
4312          """
4313          if format=="escript" or format=="str"  or format=="text":
4314             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4315          else:
4316             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4317    
4318       def substitute(self,argvals):
4319          """
4320          assigns new values to symbols in the definition of the symbol.
4321          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4322    
4323          @param argvals: new values assigned to symbols
4324          @type argvals: C{dict} with keywords of type L{Symbol}.
4325          @return: result of the substitution process. Operations are executed as much as possible.
4326          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4327          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4328          """
4329          if argvals.has_key(self):
4330             arg=argvals[self]
4331             if self.isAppropriateValue(arg):
4332                return arg
4333             else:
4334                raise TypeError,"%s: new value is not appropriate."%str(self)
4335          else:
4336             arg=self.getSubstitutedArguments(argvals)
4337             return grad(arg[0],where=arg[1])
4338    
4339       def diff(self,arg):
4340          """
4341          differential of this object
4342    
4343          @param arg: the derivative is calculated with respect to arg
4344          @type arg: L{escript.Symbol}
4345          @return: derivative with respect to C{arg}
4346          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4347          """
4348          if arg==self:
4349             return identity(self.getShape())
4350          else:
4351             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4352    
4353  def integrate(arg,where=None):  def integrate(arg,where=None):
4354      """      """
4355      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}
4356      its domain.      before integration.
4357    
4358      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4359        @type arg: L{escript.Data} or L{Symbol}
4360      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4361                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4362        @type where: C{None} or L{escript.FunctionSpace}
4363        @return: integral of arg.
4364        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4365      """      """
4366      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):  
4367         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4368      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4369         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
# Line 3658  def integrate(arg,where=None): Line 4374  def integrate(arg,where=None):
4374      else:      else:
4375        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4376    
4377    class Integrate_Symbol(DependendSymbol):
4378       """
4379       L{Symbol} representing the result of the spatial integration operator
4380       """
4381       def __init__(self,arg,where=None):
4382          """
4383          initialization of integration L{Symbol} with argument arg
4384          @param arg: argument of the integration
4385          @type arg: L{Symbol}.
4386          @param where: FunctionSpace in which the integration will be calculated.
4387                      If not present or C{None} an appropriate default is used.
4388          @type where: C{None} or L{escript.FunctionSpace}
4389          """
4390          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4391    
4392       def getMyCode(self,argstrs,format="escript"):
4393          """
4394          returns a program code that can be used to evaluate the symbol.
4395    
4396          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4397          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4398          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4399          @type format: C{str}
4400          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4401          @rtype: C{str}
4402          @raise NotImplementedError: if the requested format is not available
4403          """
4404          if format=="escript" or format=="str"  or format=="text":
4405             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4406          else:
4407             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4408    
4409       def substitute(self,argvals):
4410          """
4411          assigns new values to symbols in the definition of the symbol.
4412          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4413    
4414          @param argvals: new values assigned to symbols
4415          @type argvals: C{dict} with keywords of type L{Symbol}.
4416          @return: result of the substitution process. Operations are executed as much as possible.
4417          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4418          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4419          """
4420          if argvals.has_key(self):
4421             arg=argvals[self]
4422             if self.isAppropriateValue(arg):
4423                return arg
4424             else:
4425                raise TypeError,"%s: new value is not appropriate."%str(self)
4426          else:
4427             arg=self.getSubstitutedArguments(argvals)
4428             return integrate(arg[0],where=arg[1])
4429    
4430       def diff(self,arg):
4431          """
4432          differential of this object
4433    
4434          @param arg: the derivative is calculated with respect to arg
4435          @type arg: L{escript.Symbol}
4436          @return: derivative with respect to C{arg}
4437          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4438          """
4439          if arg==self:
4440             return identity(self.getShape())
4441          else:
4442             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4443    
4444    
4445  def interpolate(arg,where):  def interpolate(arg,where):
4446      """      """
4447      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
4448    
4449      @param arg:    interpolant      @param arg: interpolant
4450      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
4451        @param where: FunctionSpace to be interpolated to
4452        @type where: L{escript.FunctionSpace}
4453        @return: interpolated argument
4454        @rtype: C{escript.Data} or L{Symbol}
4455      """      """
4456      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4457         return Interpolated_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
4458      else:      else:
4459         return escript.Data(arg,where)         return escript.Data(arg,where)
4460    
4461  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
4462      """     """
4463      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
4464       """
4465       def __init__(self,arg,where):
4466          """
4467          initialization of interpolation L{Symbol} with argument arg
4468          @param arg: argument of the interpolation
4469          @type arg: L{Symbol}.
4470          @param where: FunctionSpace into which the argument is interpolated.
4471          @type where: L{escript.FunctionSpace}
4472          """
4473          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4474    
4475      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
4476                    be calculated.        """
4477      @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)  
4478    
4479  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
4480      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4481      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.
4482          @type format: C{str}
4483          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4484          @rtype: C{str}
4485          @raise NotImplementedError: if the requested format is not available
4486          """
4487          if format=="escript" or format=="str"  or format=="text":
4488             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
4489          else:
4490             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4491    
4492      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
4493                    to be calculated.        """
4494      """        assigns new values to symbols in the definition of the symbol.
4495      d=arg.getDomain()        The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
     return arg.interpolate(escript.FunctionOnContactOne(d))-arg.interpolate(escript.FunctionOnContactZero(d))  
4496    
4497  #=============================        @param argvals: new values assigned to symbols
4498  #        @type argvals: C{dict} with keywords of type L{Symbol}.
4499  # 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.
4500  # 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
4501  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4502          """
4503          if argvals.has_key(self):
4504             arg=argvals[self]
4505             if self.isAppropriateValue(arg):
4506                return arg
4507             else:
4508                raise TypeError,"%s: new value is not appropriate."%str(self)
4509          else:
4510             arg=self.getSubstitutedArguments(argvals)
4511             return interpolate(arg[0],where=arg[1])
4512    
4513       def diff(self,arg):
4514          """
4515          differential of this object
4516    
4517          @param arg: the derivative is calculated with respect to arg
4518          @type arg: L{escript.Symbol}
4519          @return: derivative with respect to C{arg}
4520          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
4521          """
4522          if arg==self:
4523             return identity(self.getShape())
4524          else:
4525             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4526    
 # functions involving the underlying Domain:  
4527    
4528  def transpose(arg,axis=None):  def div(arg,where=None):
4529      """      """
4530      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
4531    
4532      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
4533        @type arg: L{escript.Data} or L{Symbol}
4534        @param where: FunctionSpace in which the divergence will be calculated.
4535                      If not present or C{None} an appropriate default is used.
4536        @type where: C{None} or L{escript.FunctionSpace}
4537        @return: divergence of arg.
4538        @rtype: L{escript.Data} or L{Symbol}
4539      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
4540      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4541         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
4542      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4543         # 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)  
4544      else:      else:
4545         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
4546        if not arg.getShape()==(dim,):
4547          raise ValueError,"div: expected shape is (%s,)"%dim
4548        return trace(grad(arg,where))
4549    
4550  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
4551      """      """
4552      Return      returns the jump of arg across the continuity of the domain
4553    
4554      @param arg:      @param arg: argument
4555        @type arg: L{escript.Data} or L{Symbol}
4556        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
4557                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
4558        @type domain: C{None} or L{escript.Domain}
4559        @return: jump of arg
4560        @rtype: L{escript.Data} or L{Symbol}
4561      """      """
4562      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
4563         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)  
4564    
4565    def L2(arg):
4566        """
4567        returns the L2 norm of arg at where
4568        
4569        @param arg: function which L2 to be calculated.
4570        @type arg: L{escript.Data} or L{Symbol}
4571        @return: L2 norm of arg.
4572        @rtype: L{float} or L{Symbol}
4573        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
4574        """
4575        return sqrt(integrate(inner(arg,arg)))
4576    #=============================
4577    #
4578    
4579  def reorderComponents(arg,index):  def reorderComponents(arg,index):
4580      """      """
4581      resorts the component of arg according to index      resorts the component of arg according to index
4582    
4583      """      """
4584      pass      raise NotImplementedError
4585  #  #
4586  # $Log: util.py,v $  # $Log: util.py,v $
4587  # 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.396  
changed lines
  Added in v.775

  ViewVC Help
Powered by ViewVC 1.1.26