/[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 429 by gross, Wed Jan 11 05:53:40 2006 UTC revision 720 by gross, Thu Apr 27 10:16:05 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 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 65  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 93  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 122  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 140  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 160  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 175  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 188  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 246  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 335  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 357  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 391  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 417  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 461  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 496  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 540  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 549  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 573  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 602  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 692  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 700  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 709  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 726  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 737  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 748  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 759  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 770  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 781  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 792  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 803  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 814  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 825  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 880  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 908  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))*1.        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 954  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 990  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))*1.        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 1036  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 1072  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))*1.        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 1102  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 1134  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 1183  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 1217  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 1250  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 1288  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 1340  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 1378  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 1430  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 1468  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 1520  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 1558  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 1610  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 1648  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 1700  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 1738  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 1790  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 1828  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 1880  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 1918  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 1970  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 2008  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 2060  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 2098  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 2150  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 2188  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 2240  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 2278  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 2330  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 2368  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 2420  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 2458  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 2510  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 2548  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 2600  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 2648  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 2700  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 2741  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 2777  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 2818  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 2854  def length(arg): Line 2920  def length(arg):
2920    
2921     @param arg: argument     @param arg: argument
2922     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2923     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2924     """     """
2925     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2926    
# Line 2883  def trace(arg,axis_offset=0): Line 2949  def trace(arg,axis_offset=0):
2949        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
2950          raise ValueError,"trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"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))        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
2952        out=numarray.zeros([s1,s2],numarray.Float)        out=numarray.zeros([s1,s2],numarray.Float64)
2953        for i1 in range(s1):        for i1 in range(s1):
2954          for i2 in range(s2):          for i2 in range(s2):
2955              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
# Line 2971  class Trace_Symbol(DependendSymbol): Line 3037  class Trace_Symbol(DependendSymbol):
3037        @type format: C{str}        @type format: C{str}
3038        @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.
3039        @rtype: C{str}        @rtype: C{str}
3040        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3041        """        """
3042        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3043           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3013  class Trace_Symbol(DependendSymbol): Line 3079  class Trace_Symbol(DependendSymbol):
3079        else:        else:
3080           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3081    
3082    def transpose(arg,axis_offset=None):
3083       """
3084       returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.
3085    
3086       @param arg: argument
3087       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3088       @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.
3089                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3090       @type axis_offset: C{int}
3091       @return: transpose of arg
3092       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3093       """
3094       if isinstance(arg,numarray.NumArray):
3095          if axis_offset==None: axis_offset=int(arg.rank/2)
3096          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3097       elif isinstance(arg,escript.Data):
3098          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3099          return escript_transpose(arg,axis_offset)
3100       elif isinstance(arg,float):
3101          if not ( axis_offset==0 or axis_offset==None):
3102            raise ValueError,"transpose: axis_offset must be 0 for float argument"
3103          return arg
3104       elif isinstance(arg,int):
3105          if not ( axis_offset==0 or axis_offset==None):
3106            raise ValueError,"transpose: axis_offset must be 0 for int argument"
3107          return float(arg)
3108       elif isinstance(arg,Symbol):
3109          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3110          return Transpose_Symbol(arg,axis_offset)
3111       else:
3112          raise TypeError,"transpose: Unknown argument type."
3113    
3114    def escript_transpose(arg,axis_offset): # this should be escript._transpose
3115          "arg si a Data objects!!!"
3116          r=arg.getRank()
3117          if axis_offset<0 or axis_offset>r:
3118            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3119          s=arg.getShape()
3120          s_out=s[axis_offset:]+s[:axis_offset]
3121          out=escript.Data(0.,s_out,arg.getFunctionSpace())
3122          if r==4:
3123             if axis_offset==1:
3124                for i0 in range(s_out[0]):
3125                   for i1 in range(s_out[1]):
3126                      for i2 in range(s_out[2]):
3127                         for i3 in range(s_out[3]):
3128                             out[i0,i1,i2,i3]=arg[i3,i0,i1,i2]
3129             elif axis_offset==2:
3130                for i0 in range(s_out[0]):
3131                   for i1 in range(s_out[1]):
3132                      for i2 in range(s_out[2]):
3133                         for i3 in range(s_out[3]):
3134                             out[i0,i1,i2,i3]=arg[i2,i3,i0,i1]
3135             elif axis_offset==3:
3136                for i0 in range(s_out[0]):
3137                   for i1 in range(s_out[1]):
3138                      for i2 in range(s_out[2]):
3139                         for i3 in range(s_out[3]):
3140                             out[i0,i1,i2,i3]=arg[i1,i2,i3,i0]
3141             else:
3142                for i0 in range(s_out[0]):
3143                   for i1 in range(s_out[1]):
3144                      for i2 in range(s_out[2]):
3145                         for i3 in range(s_out[3]):
3146                             out[i0,i1,i2,i3]=arg[i0,i1,i2,i3]
3147          elif r==3:
3148             if axis_offset==1:
3149                for i0 in range(s_out[0]):
3150                   for i1 in range(s_out[1]):
3151                      for i2 in range(s_out[2]):
3152                             out[i0,i1,i2]=arg[i2,i0,i1]
3153             elif axis_offset==2:
3154                for i0 in range(s_out[0]):
3155                   for i1 in range(s_out[1]):
3156                      for i2 in range(s_out[2]):
3157                             out[i0,i1,i2]=arg[i1,i2,i0]
3158             else:
3159                for i0 in range(s_out[0]):
3160                   for i1 in range(s_out[1]):
3161                      for i2 in range(s_out[2]):
3162                             out[i0,i1,i2]=arg[i0,i1,i2]
3163          elif r==2:
3164             if axis_offset==1:
3165                for i0 in range(s_out[0]):
3166                   for i1 in range(s_out[1]):
3167                             out[i0,i1]=arg[i1,i0]
3168             else:
3169                for i0 in range(s_out[0]):
3170                   for i1 in range(s_out[1]):
3171                             out[i0,i1]=arg[i0,i1]
3172          elif r==1:
3173              for i0 in range(s_out[0]):
3174                   out[i0]=arg[i0]
3175          elif r==0:
3176                 out=arg+0.
3177          return out
3178    class Transpose_Symbol(DependendSymbol):
3179       """
3180       L{Symbol} representing the result of the transpose function
3181       """
3182       def __init__(self,arg,axis_offset=None):
3183          """
3184          initialization of transpose L{Symbol} with argument arg
3185    
3186          @param arg: argument of function
3187          @type arg: L{Symbol}.
3188          @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.
3189                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3190          @type axis_offset: C{int}
3191          """
3192          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3193          if axis_offset<0 or axis_offset>arg.getRank():
3194            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3195          s=arg.getShape()
3196          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3197    
3198       def getMyCode(self,argstrs,format="escript"):
3199          """
3200          returns a program code that can be used to evaluate the symbol.
3201    
3202          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3203          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3204          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3205          @type format: C{str}
3206          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3207          @rtype: C{str}
3208          @raise NotImplementedError: if the requested format is not available
3209          """
3210          if format=="escript" or format=="str"  or format=="text":
3211             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3212          else:
3213             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3214    
3215       def substitute(self,argvals):
3216          """
3217          assigns new values to symbols in the definition of the symbol.
3218          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3219    
3220          @param argvals: new values assigned to symbols
3221          @type argvals: C{dict} with keywords of type L{Symbol}.
3222          @return: result of the substitution process. Operations are executed as much as possible.
3223          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3224          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3225          """
3226          if argvals.has_key(self):
3227             arg=argvals[self]
3228             if self.isAppropriateValue(arg):
3229                return arg
3230             else:
3231                raise TypeError,"%s: new value is not appropriate."%str(self)
3232          else:
3233             arg=self.getSubstitutedArguments(argvals)
3234             return transpose(arg[0],axis_offset=arg[1])
3235    
3236       def diff(self,arg):
3237          """
3238          differential of this object
3239    
3240          @param arg: the derivative is calculated with respect to arg
3241          @type arg: L{escript.Symbol}
3242          @return: derivative with respect to C{arg}
3243          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3244          """
3245          if arg==self:
3246             return identity(self.getShape())
3247          else:
3248             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3249    def symmetric(arg):
3250        """
3251        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3252    
3253        @param arg: square matrix. Must have rank 2 or 4 and be square.
3254        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3255        @return: symmetric part of arg
3256        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3257        """
3258        if isinstance(arg,numarray.NumArray):
3259          if arg.rank==2:
3260            if not (arg.shape[0]==arg.shape[1]):
3261               raise ValueError,"symmetric: argument must be square."
3262          elif arg.rank==4:
3263            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3264               raise ValueError,"symmetric: argument must be square."
3265          else:
3266            raise ValueError,"symmetric: rank 2 or 4 is required."
3267          return (arg+transpose(arg))/2
3268        elif isinstance(arg,escript.Data):
3269          return escript_symmetric(arg)
3270        elif isinstance(arg,float):
3271          return arg
3272        elif isinstance(arg,int):
3273          return float(arg)
3274        elif isinstance(arg,Symbol):
3275          if arg.getRank()==2:
3276            if not (arg.getShape()[0]==arg.getShape()[1]):
3277               raise ValueError,"symmetric: argument must be square."
3278          elif arg.getRank()==4:
3279            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3280               raise ValueError,"symmetric: argument must be square."
3281          else:
3282            raise ValueError,"symmetric: rank 2 or 4 is required."
3283          return (arg+transpose(arg))/2
3284        else:
3285          raise TypeError,"symmetric: Unknown argument type."
3286    
3287    def escript_symmetric(arg): # this should be implemented in c++
3288          if arg.getRank()==2:
3289            if not (arg.getShape()[0]==arg.getShape()[1]):
3290               raise ValueError,"escript_symmetric: argument must be square."
3291            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3292            for i0 in range(arg.getShape()[0]):
3293               for i1 in range(arg.getShape()[1]):
3294                  out[i0,i1]=(arg[i0,i1]+arg[i1,i0])/2.
3295          elif arg.getRank()==4:
3296            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3297               raise ValueError,"escript_symmetric: argument must be square."
3298            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3299            for i0 in range(arg.getShape()[0]):
3300               for i1 in range(arg.getShape()[1]):
3301                  for i2 in range(arg.getShape()[2]):
3302                     for i3 in range(arg.getShape()[3]):
3303                         out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]+arg[i2,i3,i0,i1])/2.
3304          else:
3305            raise ValueError,"escript_symmetric: rank 2 or 4 is required."
3306          return out
3307    
3308    def nonsymmetric(arg):
3309        """
3310        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3311    
3312        @param arg: square matrix. Must have rank 2 or 4 and be square.
3313        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3314        @return: nonsymmetric part of arg
3315        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3316        """
3317        if isinstance(arg,numarray.NumArray):
3318          if arg.rank==2:
3319            if not (arg.shape[0]==arg.shape[1]):
3320               raise ValueError,"nonsymmetric: argument must be square."
3321          elif arg.rank==4:
3322            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3323               raise ValueError,"nonsymmetric: argument must be square."
3324          else:
3325            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3326          return (arg-transpose(arg))/2
3327        elif isinstance(arg,escript.Data):
3328          return escript_nonsymmetric(arg)
3329        elif isinstance(arg,float):
3330          return arg
3331        elif isinstance(arg,int):
3332          return float(arg)
3333        elif isinstance(arg,Symbol):
3334          if arg.getRank()==2:
3335            if not (arg.getShape()[0]==arg.getShape()[1]):
3336               raise ValueError,"nonsymmetric: argument must be square."
3337          elif arg.getRank()==4:
3338            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3339               raise ValueError,"nonsymmetric: argument must be square."
3340          else:
3341            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3342          return (arg-transpose(arg))/2
3343        else:
3344          raise TypeError,"nonsymmetric: Unknown argument type."
3345    
3346    def escript_nonsymmetric(arg): # this should be implemented in c++
3347          if arg.getRank()==2:
3348            if not (arg.getShape()[0]==arg.getShape()[1]):
3349               raise ValueError,"escript_nonsymmetric: argument must be square."
3350            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3351            for i0 in range(arg.getShape()[0]):
3352               for i1 in range(arg.getShape()[1]):
3353                  out[i0,i1]=(arg[i0,i1]-arg[i1,i0])/2.
3354          elif arg.getRank()==4:
3355            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3356               raise ValueError,"escript_nonsymmetric: argument must be square."
3357            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3358            for i0 in range(arg.getShape()[0]):
3359               for i1 in range(arg.getShape()[1]):
3360                  for i2 in range(arg.getShape()[2]):
3361                     for i3 in range(arg.getShape()[3]):
3362                         out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]-arg[i2,i3,i0,i1])/2.
3363          else:
3364            raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."
3365          return out
3366    
3367    
3368    def inverse(arg):
3369        """
3370        returns the inverse of the square matrix arg.
3371    
3372        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3373        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3374        @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3375        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3376        @note: for L{escript.Data} objects the dimension is restricted to 3.
3377        """
3378        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3379        if isinstance(arg,numarray.NumArray):
3380          return numarray.linear_algebra.inverse(arg)
3381        elif isinstance(arg,escript.Data):
3382          return escript_inverse(arg)
3383        elif isinstance(arg,float):
3384          return 1./arg
3385        elif isinstance(arg,int):
3386          return 1./float(arg)
3387        elif isinstance(arg,Symbol):
3388          return Inverse_Symbol(arg)
3389        else:
3390          raise TypeError,"inverse: Unknown argument type."
3391    
3392    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3393          "arg is a Data objects!!!"
3394          if not arg.getRank()==2:
3395            raise ValueError,"escript_inverse: argument must have rank 2"
3396          s=arg.getShape()      
3397          if not s[0] == s[1]:
3398            raise ValueError,"escript_inverse: argument must be a square matrix."
3399          out=escript.Data(0.,s,arg.getFunctionSpace())
3400          if s[0]==1:
3401              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3402                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3403              out[0,0]=1./arg[0,0]
3404          elif s[0]==2:
3405              A11=arg[0,0]
3406              A12=arg[0,1]
3407              A21=arg[1,0]
3408              A22=arg[1,1]
3409              D = A11*A22-A12*A21
3410              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3411                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3412              D=1./D
3413              out[0,0]= A22*D
3414              out[1,0]=-A21*D
3415              out[0,1]=-A12*D
3416              out[1,1]= A11*D
3417          elif s[0]==3:
3418              A11=arg[0,0]
3419              A21=arg[1,0]
3420              A31=arg[2,0]
3421              A12=arg[0,1]
3422              A22=arg[1,1]
3423              A32=arg[2,1]
3424              A13=arg[0,2]
3425              A23=arg[1,2]
3426              A33=arg[2,2]
3427              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3428              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3429                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3430              D=1./D
3431              out[0,0]=(A22*A33-A23*A32)*D
3432              out[1,0]=(A31*A23-A21*A33)*D
3433              out[2,0]=(A21*A32-A31*A22)*D
3434              out[0,1]=(A13*A32-A12*A33)*D
3435              out[1,1]=(A11*A33-A31*A13)*D
3436              out[2,1]=(A12*A31-A11*A32)*D
3437              out[0,2]=(A12*A23-A13*A22)*D
3438              out[1,2]=(A13*A21-A11*A23)*D
3439              out[2,2]=(A11*A22-A12*A21)*D
3440          else:
3441             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3442          return out
3443    
3444    class Inverse_Symbol(DependendSymbol):
3445       """
3446       L{Symbol} representing the result of the inverse function
3447       """
3448       def __init__(self,arg):
3449          """
3450          initialization of inverse L{Symbol} with argument arg
3451          @param arg: argument of function
3452          @type arg: L{Symbol}.
3453          """
3454          if not arg.getRank()==2:
3455            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3456          s=arg.getShape()
3457          if not s[0] == s[1]:
3458            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3459          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3460    
3461       def getMyCode(self,argstrs,format="escript"):
3462          """
3463          returns a program code that can be used to evaluate the symbol.
3464    
3465          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3466          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3467          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3468          @type format: C{str}
3469          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3470          @rtype: C{str}
3471          @raise NotImplementedError: if the requested format is not available
3472          """
3473          if format=="escript" or format=="str"  or format=="text":
3474             return "inverse(%s)"%argstrs[0]
3475          else:
3476             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3477    
3478       def substitute(self,argvals):
3479          """
3480          assigns new values to symbols in the definition of the symbol.
3481          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3482    
3483          @param argvals: new values assigned to symbols
3484          @type argvals: C{dict} with keywords of type L{Symbol}.
3485          @return: result of the substitution process. Operations are executed as much as possible.
3486          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3487          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3488          """
3489          if argvals.has_key(self):
3490             arg=argvals[self]
3491             if self.isAppropriateValue(arg):
3492                return arg
3493             else:
3494                raise TypeError,"%s: new value is not appropriate."%str(self)
3495          else:
3496             arg=self.getSubstitutedArguments(argvals)
3497             return inverse(arg[0])
3498    
3499       def diff(self,arg):
3500          """
3501          differential of this object
3502    
3503          @param arg: the derivative is calculated with respect to arg
3504          @type arg: L{escript.Symbol}
3505          @return: derivative with respect to C{arg}
3506          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3507          """
3508          if arg==self:
3509             return identity(self.getShape())
3510          else:
3511             return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)
3512    
3513    def eigenvalues(arg):
3514        """
3515        returns the eigenvalues of the square matrix arg.
3516    
3517        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3518                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3519        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3520        @return: the eigenvalues in increasing order.
3521        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3522        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3523        """
3524        if isinstance(arg,numarray.NumArray):
3525          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3526          out.sort()
3527          return out
3528        elif isinstance(arg,escript.Data):
3529          return arg._eigenvalues()
3530        elif isinstance(arg,Symbol):
3531          if not arg.getRank()==2:
3532            raise ValueError,"eigenvalues: argument must have rank 2"
3533          s=arg.getShape()      
3534          if not s[0] == s[1]:
3535            raise ValueError,"eigenvalues: argument must be a square matrix."
3536          if s[0]==1:
3537              return arg[0]
3538          elif s[0]==2:
3539              arg1=symmetric(arg)
3540              A11=arg1[0,0]
3541              A12=arg1[0,1]
3542              A22=arg1[1,1]
3543              trA=(A11+A22)/2.
3544              A11-=trA
3545              A22-=trA
3546              s=sqrt(A12**2-A11*A22)
3547              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3548          elif s[0]==3:
3549              arg1=symmetric(arg)
3550              A11=arg1[0,0]
3551              A12=arg1[0,1]
3552              A22=arg1[1,1]
3553              A13=arg1[0,2]
3554              A23=arg1[1,2]
3555              A33=arg1[2,2]
3556              trA=(A11+A22+A33)/3.
3557              A11-=trA
3558              A22-=trA
3559              A33-=trA
3560              A13_2=A13**2
3561              A23_2=A23**2
3562              A12_2=A12**2
3563              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3564              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3565              sq_p=sqrt(p/3.)
3566              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
3567              sq_p*=2.
3568              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3569               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3570               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3571              return trA+sq_p*f
3572          else:
3573             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3574        elif isinstance(arg,float):
3575          return arg
3576        elif isinstance(arg,int):
3577          return float(arg)
3578        else:
3579          raise TypeError,"eigenvalues: Unknown argument type."
3580    
3581    def eigenvalues_and_eigenvectors(arg):
3582        """
3583        returns the eigenvalues and eigenvectors of the square matrix arg.
3584    
3585        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3586                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3587        @type arg: L{escript.Data}
3588        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3589                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3590                 the eigenvector coresponding to the i-th eigenvalue.
3591        @rtype: L{tuple} of L{escript.Data}.
3592        @note: The dimension is restricted to 3.
3593        """
3594        if isinstance(arg,numarray.NumArray):
3595          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3596        elif isinstance(arg,escript.Data):
3597          return arg._eigenvalues_and_eigenvectors()
3598        elif isinstance(arg,Symbol):
3599          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3600        elif isinstance(arg,float):
3601          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3602        elif isinstance(arg,int):
3603          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3604        else:
3605          raise TypeError,"eigenvalues: Unknown argument type."
3606  #=======================================================  #=======================================================
3607  #  Binary operations:  #  Binary operations:
3608  #=======================================================  #=======================================================
# Line 3072  class Add_Symbol(DependendSymbol): Line 3662  class Add_Symbol(DependendSymbol):
3662        @type format: C{str}        @type format: C{str}
3663        @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.
3664        @rtype: C{str}        @rtype: C{str}
3665        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3666        """        """
3667        if format=="str" or format=="text":        if format=="str" or format=="text":
3668           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 3131  def mult(arg0,arg1): Line 3721  def mult(arg0,arg1):
3721         """         """
3722         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3723         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3724            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3725         else:         else:
3726            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3727                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3171  class Mult_Symbol(DependendSymbol): Line 3761  class Mult_Symbol(DependendSymbol):
3761        @type format: C{str}        @type format: C{str}
3762        @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.
3763        @rtype: C{str}        @rtype: C{str}
3764        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3765        """        """
3766        if format=="str" or format=="text":        if format=="str" or format=="text":
3767           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3231  def quotient(arg0,arg1): Line 3821  def quotient(arg0,arg1):
3821         """         """
3822         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3823         if testForZero(args[0]):         if testForZero(args[0]):
3824            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3825         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3826            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3827               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3276  class Quotient_Symbol(DependendSymbol): Line 3866  class Quotient_Symbol(DependendSymbol):
3866        @type format: C{str}        @type format: C{str}
3867        @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.
3868        @rtype: C{str}        @rtype: C{str}
3869        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3870        """        """
3871        if format=="str" or format=="text":        if format=="str" or format=="text":
3872           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3337  def power(arg0,arg1): Line 3927  def power(arg0,arg1):
3927         """         """
3928         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3929         if testForZero(args[0]):         if testForZero(args[0]):
3930            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3931         elif testForZero(args[1]):         elif testForZero(args[1]):
3932            return numarray.ones(args[0],numarray.Float)            return numarray.ones(pokeShape(args[1]),numarray.Float64)
3933         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3934            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3935         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 3380  class Power_Symbol(DependendSymbol): Line 3970  class Power_Symbol(DependendSymbol):
3970        @type format: C{str}        @type format: C{str}
3971        @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.
3972        @rtype: C{str}        @rtype: C{str}
3973        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3974        """        """
3975        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
3976           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3469  def clip(arg,minval=0.,maxval=1.): Line 4059  def clip(arg,minval=0.,maxval=1.):
4059      @param arg: argument      @param arg: argument
4060      @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}
4061      @param minval: lower range      @param minval: lower range
4062      @type arg: C{float}      @type minval: C{float}
4063      @param maxval: upper range      @param maxval: upper range
4064      @type arg: C{float}      @type maxval: C{float}
4065      @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
4066               less then maxval are unchanged.               less then maxval are unchanged.
4067      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
# Line 3496  def inner(arg0,arg1): Line 4086  def inner(arg0,arg1):
4086      @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}
4087      @param arg1: second argument      @param arg1: second argument
4088      @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}
4089      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4090      @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
4091      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4092      """      """
# Line 3512  def matrixmult(arg0,arg1): Line 4102  def matrixmult(arg0,arg1):
4102    
4103      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4104    
4105            or      or
4106    
4107      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4108    
# Line 3540  def outer(arg0,arg1): Line 4130  def outer(arg0,arg1):
4130    
4131      out[t,s]=arg0[t]*arg1[s]      out[t,s]=arg0[t]*arg1[s]
4132    
4133      where s runs through arg0.Shape      where
4134            t runs through arg1.Shape  
4135            - s runs through arg0.Shape
4136            - t runs through arg1.Shape
4137    
4138      @param arg0: first argument      @param arg0: first argument
4139      @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 3556  def outer(arg0,arg1): Line 4148  def outer(arg0,arg1):
4148  def tensormult(arg0,arg1):  def tensormult(arg0,arg1):
4149      """      """
4150      the tensor product of the two argument:      the tensor product of the two argument:
   
4151            
4152      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4153    
4154      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4155    
4156                   or      or
4157    
4158      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4159    
# Line 3571  def tensormult(arg0,arg1): Line 4162  def tensormult(arg0,arg1):
4162    
4163      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]
4164                                
4165                   or      or
4166    
4167      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]
4168    
4169                   or      or
4170    
4171      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]
4172    
# Line 3604  def generalTensorProduct(arg0,arg1,axis_ Line 4195  def generalTensorProduct(arg0,arg1,axis_
4195    
4196      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4197    
4198      where s runs through arg0.Shape[:arg0.Rank-axis_offset]      where
4199            r runs trough arg0.Shape[:axis_offset]  
4200            t runs through arg1.Shape[axis_offset:]          - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4201            - r runs trough arg0.Shape[:axis_offset]
4202            - t runs through arg1.Shape[axis_offset:]
4203    
4204      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  
4205      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 3636  def generalTensorProduct(arg0,arg1,axis_ Line 4229  def generalTensorProduct(arg0,arg1,axis_
4229             for i in sh1[:axis_offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4230             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4231             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4232             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4233             for i0 in range(d0):             for i0 in range(d0):
4234                      for i1 in range(d1):                      for i1 in range(d1):
4235                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
# Line 3685  class GeneralTensorProduct_Symbol(Depend Line 4278  class GeneralTensorProduct_Symbol(Depend
4278        @type format: C{str}        @type format: C{str}
4279        @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.
4280        @rtype: C{str}        @rtype: C{str}
4281        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4282        """        """
4283        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4284           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
# Line 3776  def grad(arg,where=None): Line 4369  def grad(arg,where=None):
4369                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4370      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4371      @return: gradient of arg.      @return: gradient of arg.
4372      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4373      """      """
4374      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4375         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3804  class Grad_Symbol(DependendSymbol): Line 4397  class Grad_Symbol(DependendSymbol):
4397        d=arg.getDim()        d=arg.getDim()
4398        if d==None:        if d==None:
4399           raise ValueError,"argument must have a spatial dimension"           raise ValueError,"argument must have a spatial dimension"
4400        super(Grad_Symbol,self).__init__(args=[arg,where],shape=tuple(list(arg.getShape()).extend(d)),dim=d)        super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4401    
4402     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4403        """        """
# Line 3816  class Grad_Symbol(DependendSymbol): Line 4409  class Grad_Symbol(DependendSymbol):
4409        @type format: C{str}        @type format: C{str}
4410        @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.
4411        @rtype: C{str}        @rtype: C{str}
4412        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4413        """        """
4414        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4415           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 3869  def integrate(arg,where=None): Line 4462  def integrate(arg,where=None):
4462                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4463      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4464      @return: integral of arg.      @return: integral of arg.
4465      @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}      @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4466      """      """
4467      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4468         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
# Line 3907  class Integrate_Symbol(DependendSymbol): Line 4500  class Integrate_Symbol(DependendSymbol):
4500        @type format: C{str}        @type format: C{str}
4501        @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.
4502        @rtype: C{str}        @rtype: C{str}
4503        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4504        """        """
4505        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4506           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 3959  def interpolate(arg,where): Line 4552  def interpolate(arg,where):
4552      @param where: FunctionSpace to be interpolated to      @param where: FunctionSpace to be interpolated to
4553      @type where: L{escript.FunctionSpace}      @type where: L{escript.FunctionSpace}
4554      @return: interpolated argument      @return: interpolated argument
4555      @rtype:  C{escript.Data} or L{Symbol}      @rtype: C{escript.Data} or L{Symbol}
4556      """      """
4557      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4558         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
# Line 3990  class Interpolate_Symbol(DependendSymbol Line 4583  class Interpolate_Symbol(DependendSymbol
4583        @type format: C{str}        @type format: C{str}
4584        @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.
4585        @rtype: C{str}        @rtype: C{str}
4586        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4587        """        """
4588        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4589           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4043  def div(arg,where=None): Line 4636  def div(arg,where=None):
4636                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4637      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4638      @return: divergence of arg.      @return: divergence of arg.
4639      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4640      """      """
4641      if not arg.getShape()==(arg.getDim(),):      if isinstance(arg,Symbol):
4642        raise ValueError,"div: expected shape is (%s,)"%arg.getDim()          dim=arg.getDim()
4643        elif isinstance(arg,escript.Data):
4644            dim=arg.getDomain().getDim()
4645        else:
4646            raise TypeError,"div: argument type not supported"
4647        if not arg.getShape()==(dim,):
4648          raise ValueError,"div: expected shape is (%s,)"%dim
4649      return trace(grad(arg,where))      return trace(grad(arg,where))
4650    
4651  def jump(arg,domain=None):  def jump(arg,domain=None):
# Line 4059  def jump(arg,domain=None): Line 4658  def jump(arg,domain=None):
4658                     the domain of arg is used. If arg is a L{Symbol} the domain must be present.                     the domain of arg is used. If arg is a L{Symbol} the domain must be present.
4659      @type domain: C{None} or L{escript.Domain}      @type domain: C{None} or L{escript.Domain}
4660      @return: jump of arg      @return: jump of arg
4661      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4662      """      """
4663      if domain==None: domain=arg.getDomain()      if domain==None: domain=arg.getDomain()
4664      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
 #=============================  
 #  
 # wrapper for various functions: if the argument has attribute the function name  
 # as an argument it calls the corresponding methods. Otherwise the corresponding  
 # numarray function is called.  
   
 # functions involving the underlying Domain:  
4665    
4666  def transpose(arg,axis=None):  def L2(arg):
4667      """      """
4668      Returns the transpose of the Data object arg.      returns the L2 norm of arg at where
4669        
4670      @param arg:      @param arg: function which L2 to be calculated.
4671        @type arg: L{escript.Data} or L{Symbol}
4672        @return: L2 norm of arg.
4673        @rtype: L{float} or L{Symbol}
4674        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
4675      """      """
4676      if axis==None:      return sqrt(integrate(inner(arg,arg)))
4677         r=0  #=============================
4678         if hasattr(arg,"getRank"): r=arg.getRank()  #
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
     if isinstance(arg,Symbol):  
        return Transpose_Symbol(arg,axis=r)  
     if isinstance(arg,escript.Data):  
        # hack for transpose  
        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)  
     else:  
        return numarray.transpose(arg,axis=axis)  
   
   
4679    
4680  def reorderComponents(arg,index):  def reorderComponents(arg,index):
4681      """      """
4682      resorts the component of arg according to index      resorts the component of arg according to index
4683    
4684      """      """
4685      pass      raise NotImplementedError
4686  #  #
4687  # $Log: util.py,v $  # $Log: util.py,v $
4688  # 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.429  
changed lines
  Added in v.720

  ViewVC Help
Powered by ViewVC 1.1.26