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

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

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

revision 341 by gross, Mon Dec 12 05:26:10 2005 UTC revision 881 by gross, Thu Oct 26 02:54:47 2006 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2    
3  """  """
4  Utility functions for escript  Utility functions for escript
5    
 @remark:  This module is under construction and is still tested!!!  
   
6  @var __author__: name of author  @var __author__: name of author
7  @var __licence__: licence agreement  @var __copyright__: copyrights
8    @var __license__: licence agreement
9  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
10  @var __version__: version  @var __version__: version
11  @var __date__: date of the version  @var __date__: date of the version
12  """  """
13                                                                                                                                                                                                                                                                                                                                                                                                            
14  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
15  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                        http://www.access.edu.au
17                    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19                 http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
21  __version__="$Revision: 329 $"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
23    
24    
# Line 32  import math Line 26  import math
26  import numarray  import numarray
27  import escript  import escript
28  import os  import os
29    from esys.escript import C_GeneralTensorProduct
 # missing tests:  
   
 # def pokeShape(arg):  
 # def pokeDim(arg):  
 # def commonShape(arg0,arg1):  
 # def commonDim(*args):  
 # def testForZero(arg):  
 # def matchType(arg0=0.,arg1=0.):  
 # def matchShape(arg0,arg1):  
   
 # def maximum(arg0,arg1):  
 # def minimum(arg0,arg1):  
   
 # def transpose(arg,axis=None):  
 # def trace(arg,axis0=0,axis1=1):  
 # def reorderComponents(arg,index):  
   
 # def integrate(arg,where=None):  
 # def interpolate(arg,where):  
 # def div(arg,where=None):  
 # def grad(arg,where=None):  
   
 #  
 # slicing: get  
 #          set  
 #  
 # and derivatives  
30    
31  #=========================================================  #=========================================================
32  #   some helpers:  #   some helpers:
# Line 68  def saveVTK(filename,domain=None,**data) Line 35  def saveVTK(filename,domain=None,**data)
35      """      """
36      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.
37    
38      Example:      Example::
39    
40         tmp=Scalar(..)         tmp=Scalar(..)
41         v=Vector(..)         v=Vector(..)
# Line 96  def saveDX(filename,domain=None,**data): Line 63  def saveDX(filename,domain=None,**data):
63      """      """
64      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.
65    
66      Example:      Example::
67    
68         tmp=Scalar(..)         tmp=Scalar(..)
69         v=Vector(..)         v=Vector(..)
# Line 125  def kronecker(d=3): Line 92  def kronecker(d=3):
92     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
93    
94     @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
95     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
96     @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
97     @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}  
98     """     """
99     return identityTensor(d)     return identityTensor(d)
100    
# Line 143  def identity(shape=()): Line 109  def identity(shape=()):
109     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
110     """     """
111     if len(shape)>0:     if len(shape)>0:
112        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
113        if len(shape)==1:        if len(shape)==1:
114            for i0 in range(shape[0]):            for i0 in range(shape[0]):
115               out[i0,i0]=1.               out[i0,i0]=1.
   
116        elif len(shape)==2:        elif len(shape)==2:
117            for i0 in range(shape[0]):            for i0 in range(shape[0]):
118               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 163  def identityTensor(d=3): Line 128  def identityTensor(d=3):
128     return the dxd identity matrix     return the dxd identity matrix
129    
130     @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
131     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
132     @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
133     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
134     """     """
135     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
136        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
137     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
138           return identity((d.getDim(),))
139       else:
140           return identity((d,))
141    
142  def identityTensor4(d=3):  def identityTensor4(d=3):
143     """     """
# Line 178  def identityTensor4(d=3): Line 146  def identityTensor4(d=3):
146     @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
147     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
148     @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
149     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
150     """     """
151     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
152        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
153     return identity((d,d))     elif isinstance(d,escript.Domain):
154           return identity((d.getDim(),d.getDim()))
155       else:
156           return identity((d,d))
157    
158  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
159     """     """
# Line 191  def unitVector(i=0,d=3): Line 162  def unitVector(i=0,d=3):
162     @param i: index     @param i: index
163     @type i: C{int}     @type i: C{int}
164     @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
165     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
166     @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
167     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
168     """     """
169     return kronecker(d)[i]     return kronecker(d)[i]
170    
# Line 249  def inf(arg): Line 220  def inf(arg):
220    
221      @param arg: argument      @param arg: argument
222      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
223      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
224      @rtype: C{float}      @rtype: C{float}
225      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
226      """      """
# Line 338  def commonDim(*args): Line 309  def commonDim(*args):
309      """      """
310      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.
311    
312      @param *args: given objects      @param args: given objects
313      @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
314               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
315      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 360  def testForZero(arg): Line 331  def testForZero(arg):
331    
332      @param arg: a given object      @param arg: a given object
333      @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}
334      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
335      @rtype : C{bool}      @rtype: C{bool}
336      """      """
337      try:      if isinstance(arg,numarray.NumArray):
338           return not Lsup(arg)>0.
339        elif isinstance(arg,escript.Data):
340           return False
341        elif isinstance(arg,float):
342           return not Lsup(arg)>0.
343        elif isinstance(arg,int):
344         return not Lsup(arg)>0.         return not Lsup(arg)>0.
345      except TypeError:      elif isinstance(arg,Symbol):
346           return False
347        else:
348         return False         return False
349    
350  def matchType(arg0=0.,arg1=0.):  def matchType(arg0=0.,arg1=0.):
# Line 386  def matchType(arg0=0.,arg1=0.): Line 365  def matchType(arg0=0.,arg1=0.):
365         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
366            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
367         elif isinstance(arg1,float):         elif isinstance(arg1,float):
368            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
369         elif isinstance(arg1,int):         elif isinstance(arg1,int):
370            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
371         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
372            pass            pass
373         else:         else:
# Line 412  def matchType(arg0=0.,arg1=0.): Line 391  def matchType(arg0=0.,arg1=0.):
391         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
392            pass            pass
393         elif isinstance(arg1,float):         elif isinstance(arg1,float):
394            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
395         elif isinstance(arg1,int):         elif isinstance(arg1,int):
396            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
397         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
398            pass            pass
399         else:         else:
400            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
401      elif isinstance(arg0,float):      elif isinstance(arg0,float):
402         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
403            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
404         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
405            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
406         elif isinstance(arg1,float):         elif isinstance(arg1,float):
407            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
408            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
409         elif isinstance(arg1,int):         elif isinstance(arg1,int):
410            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
411            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
412         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
413            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
414         else:         else:
415            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
416      elif isinstance(arg0,int):      elif isinstance(arg0,int):
417         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
418            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
419         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
420            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
421         elif isinstance(arg1,float):         elif isinstance(arg1,float):
422            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
423            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
424         elif isinstance(arg1,int):         elif isinstance(arg1,int):
425            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
426            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
427         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
428            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
429         else:         else:
430            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
431      else:      else:
# Line 456  def matchType(arg0=0.,arg1=0.): Line 435  def matchType(arg0=0.,arg1=0.):
435    
436  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
437      """      """
438            return representations of arg0 amd arg1 which ahve the same shape
   
     If shape is not given the shape "largest" shape of args is used.  
439    
440      @param args: a given ob      @param arg0: a given object
441      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
442      @return: True if the argument is identical to zero.      @param arg1: a given object
443      @rtype: C{list} of C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
444        @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
445        @rtype: C{tuple}
446      """      """
447      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
448      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
449      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
450      if len(sh0)<len(sh):      if len(sh0)<len(sh):
451         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
452      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
453         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
454      else:      else:
455         return arg0,arg1         return arg0,arg1
456  #=========================================================  #=========================================================
# Line 491  class Symbol(object): Line 470  class Symbol(object):
470         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
471         symbols or any other object.         symbols or any other object.
472    
473         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
474         @type arg: C{list}         @type args: C{list}
475         @param shape: the shape         @param shape: the shape
476         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
477         @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
# Line 535  class Symbol(object): Line 514  class Symbol(object):
514         """         """
515         the shape of the symbol.         the shape of the symbol.
516    
517         @return : the shape of the symbol.         @return: the shape of the symbol.
518         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
519         """         """
520         return self.__shape         return self.__shape
# Line 544  class Symbol(object): Line 523  class Symbol(object):
523         """         """
524         the spatial dimension         the spatial dimension
525    
526         @return : the spatial dimension         @return: the spatial dimension
527         @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.
528         """         """
529         return self.__dim         return self.__dim
# Line 568  class Symbol(object): Line 547  class Symbol(object):
547         """         """
548         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.
549    
550         @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].
551         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
552         @rtype: C{list} of objects         @rtype: C{list} of objects
553         @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.
554         """         """
555         out=[]         out=[]
556         for a in self.getArgument():         for a in self.getArgument():
# Line 597  class Symbol(object): Line 576  class Symbol(object):
576            else:            else:
577                s=pokeShape(s)+arg.getShape()                s=pokeShape(s)+arg.getShape()
578                if len(s)>0:                if len(s)>0:
579                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
580                else:                else:
581                   out.append(a)                   out.append(a)
582         return out         return out
# Line 687  class Symbol(object): Line 666  class Symbol(object):
666         else:         else:
667            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
668            if len(s)>0:            if len(s)>0:
669               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
670            else:            else:
671               return 0.               return 0.
672    
# Line 695  class Symbol(object): Line 674  class Symbol(object):
674         """         """
675         returns -self.         returns -self.
676    
677         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
678         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
679         """         """
680         return self*(-1.)         return self*(-1.)
# Line 704  class Symbol(object): Line 683  class Symbol(object):
683         """         """
684         returns +self.         returns +self.
685    
686         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
687         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
688         """         """
689         return self*(1.)         return self*(1.)
690    
691     def __abs__(self):     def __abs__(self):
692         """         """
693         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
694         """         """
695         return Abs_Symbol(self)         return Abs_Symbol(self)
696    
# Line 721  class Symbol(object): Line 700  class Symbol(object):
700    
701         @param other: object to be added to this object         @param other: object to be added to this object
702         @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}.
703         @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}
704         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
705         """         """
706         return add(self,other)         return add(self,other)
# Line 732  class Symbol(object): Line 711  class Symbol(object):
711    
712         @param other: object this object is added to         @param other: object this object is added to
713         @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}.
714         @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
715         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
716         """         """
717         return add(other,self)         return add(other,self)
# Line 743  class Symbol(object): Line 722  class Symbol(object):
722    
723         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
724         @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}.
725         @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
726         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
727         """         """
728         return add(self,-other)         return add(self,-other)
# Line 754  class Symbol(object): Line 733  class Symbol(object):
733    
734         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
735         @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}.
736         @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}.
737         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
738         """         """
739         return add(-self,other)         return add(-self,other)
# Line 765  class Symbol(object): Line 744  class Symbol(object):
744    
745         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
746         @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}.
747         @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}.
748         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
749         """         """
750         return mult(self,other)         return mult(self,other)
# Line 776  class Symbol(object): Line 755  class Symbol(object):
755    
756         @param other: object this object is multiplied with         @param other: object this object is multiplied with
757         @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}.
758         @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.
759         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
760         """         """
761         return mult(other,self)         return mult(other,self)
# Line 787  class Symbol(object): Line 766  class Symbol(object):
766    
767         @param other: object dividing this object         @param other: object dividing this object
768         @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}.
769         @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}
770         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
771         """         """
772         return quotient(self,other)         return quotient(self,other)
# Line 798  class Symbol(object): Line 777  class Symbol(object):
777    
778         @param other: object dividing this object         @param other: object dividing this object
779         @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}.
780         @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
781         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
782         """         """
783         return quotient(other,self)         return quotient(other,self)
# Line 809  class Symbol(object): Line 788  class Symbol(object):
788    
789         @param other: exponent         @param other: exponent
790         @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}.
791         @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}
792         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
793         """         """
794         return power(self,other)         return power(self,other)
# Line 820  class Symbol(object): Line 799  class Symbol(object):
799    
800         @param other: basis         @param other: basis
801         @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}.
802         @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
803         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
804         """         """
805         return power(other,self)         return power(other,self)
806    
807       def __getitem__(self,index):
808           """
809           returns the slice defined by index
810    
811           @param index: defines a
812           @type index: C{slice} or C{int} or a C{tuple} of them
813           @return: a L{Symbol} representing the slice defined by index
814           @rtype: L{DependendSymbol}
815           """
816           return GetSlice_Symbol(self,index)
817    
818  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
819     """     """
820     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.
821     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  
822        
823     Example:     Example::
824        
825     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
826     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
827     print u1==u2       print u1==u2
828     False       False
829        
830        but       but::
831    
832     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
833     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
834     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
835     print u1==u2, u1==u3       print u1==u2, u1==u3
836     True False       True False
837    
838     @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.
839     """     """
# Line 875  class DependendSymbol(Symbol): Line 865  class DependendSymbol(Symbol):
865  #=========================================================  #=========================================================
866  #  Unary operations prserving the shape  #  Unary operations prserving the shape
867  #========================================================  #========================================================
868    class GetSlice_Symbol(DependendSymbol):
869       """
870       L{Symbol} representing getting a slice for a L{Symbol}
871       """
872       def __init__(self,arg,index):
873          """
874          initialization of wherePositive L{Symbol} with argument arg
875          @param arg: argument
876          @type arg: L{Symbol}.
877          @param index: defines index
878          @type index: C{slice} or C{int} or a C{tuple} of them
879          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
880          @raises ValueError: if a step is given
881          """
882          if not isinstance(index,tuple): index=(index,)
883          if len(index)>arg.getRank():
884               raise IndexError,"GetSlice_Symbol: index out of range."
885          sh=()
886          index2=()
887          for i in range(len(index)):
888             ix=index[i]
889             if isinstance(ix,int):
890                if ix<0 or ix>=arg.getShape()[i]:
891                   raise ValueError,"GetSlice_Symbol: index out of range."
892                index2=index2+(ix,)
893             else:
894               if not ix.step==None:
895                 raise ValueError,"GetSlice_Symbol: steping is not supported."
896               if ix.start==None:
897                  s=0
898               else:
899                  s=ix.start
900               if ix.stop==None:
901                  e=arg.getShape()[i]
902               else:
903                  e=ix.stop
904                  if e>arg.getShape()[i]:
905                     raise IndexError,"GetSlice_Symbol: index out of range."
906               index2=index2+(slice(s,e),)
907               if e>s:
908                   sh=sh+(e-s,)
909               elif s>e:
910                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
911          for i in range(len(index),arg.getRank()):
912              index2=index2+(slice(0,arg.getShape()[i]),)
913              sh=sh+(arg.getShape()[i],)
914          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
915    
916       def getMyCode(self,argstrs,format="escript"):
917          """
918          returns a program code that can be used to evaluate the symbol.
919    
920          @param argstrs: gives for each argument a string representing the argument for the evaluation.
921          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
922          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
923          @type format: C{str}
924          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
925          @rtype: C{str}
926          @raise NotImplementedError: if the requested format is not available
927          """
928          if format=="escript" or format=="str"  or format=="text":
929             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
930          else:
931             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
932    
933       def substitute(self,argvals):
934          """
935          assigns new values to symbols in the definition of the symbol.
936          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
937    
938          @param argvals: new values assigned to symbols
939          @type argvals: C{dict} with keywords of type L{Symbol}.
940          @return: result of the substitution process. Operations are executed as much as possible.
941          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
942          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
943          """
944          if argvals.has_key(self):
945             arg=argvals[self]
946             if self.isAppropriateValue(arg):
947                return arg
948             else:
949                raise TypeError,"%s: new value is not appropriate."%str(self)
950          else:
951             args=self.getSubstitutedArguments(argvals)
952             arg=args[0]
953             index=args[1]
954             return arg.__getitem__(index)
955    
956  def log10(arg):  def log10(arg):
957     """     """
958     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
959    
960     @param arg: argument     @param arg: argument
961     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
962     @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.
963     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
964     """     """
965     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 903  def wherePositive(arg): Line 981  def wherePositive(arg):
981    
982     @param arg: argument     @param arg: argument
983     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
984     @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.
985     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
986     """     """
987     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
988        if arg.rank==0:        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
989           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
990             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))  
991     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
992        return arg._wherePositive()        return arg._wherePositive()
993     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 953  class WherePositive_Symbol(DependendSymb Line 1027  class WherePositive_Symbol(DependendSymb
1027        @type format: C{str}        @type format: C{str}
1028        @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.
1029        @rtype: C{str}        @rtype: C{str}
1030        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1031        """        """
1032        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1033            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 989  def whereNegative(arg): Line 1063  def whereNegative(arg):
1063    
1064     @param arg: argument     @param arg: argument
1065     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1066     @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.
1067     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1068     """     """
1069     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1070        if arg.rank==0:        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1071           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1072             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))  
1073     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1074        return arg._whereNegative()        return arg._whereNegative()
1075     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1039  class WhereNegative_Symbol(DependendSymb Line 1109  class WhereNegative_Symbol(DependendSymb
1109        @type format: C{str}        @type format: C{str}
1110        @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.
1111        @rtype: C{str}        @rtype: C{str}
1112        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1113        """        """
1114        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1115            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1075  def whereNonNegative(arg): Line 1145  def whereNonNegative(arg):
1145    
1146     @param arg: argument     @param arg: argument
1147     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1148     @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.
1149     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1150     """     """
1151     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1152        if arg.rank==0:        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1153           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1154             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))  
1155     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1156        return arg._whereNonNegative()        return arg._whereNonNegative()
1157     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1109  def whereNonPositive(arg): Line 1175  def whereNonPositive(arg):
1175    
1176     @param arg: argument     @param arg: argument
1177     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1178     @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.
1179     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1180     """     """
1181     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1182        if arg.rank==0:        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1183           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1184             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.  
1185     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1186        return arg._whereNonPositive()        return arg._whereNonPositive()
1187     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1145  def whereZero(arg,tol=0.): Line 1207  def whereZero(arg,tol=0.):
1207     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1208     @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.
1209     @type tol: C{float}     @type tol: C{float}
1210     @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.
1211     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1212     """     """
1213     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1214        if arg.rank==0:        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1215           if abs(arg)<=tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1216             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
1217     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1218        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1219     elif isinstance(arg,float):     elif isinstance(arg,float):
1220        if abs(arg)<=tol:        if abs(arg)<=tol:
1221          return 1.          return 1.
# Line 1198  class WhereZero_Symbol(DependendSymbol): Line 1253  class WhereZero_Symbol(DependendSymbol):
1253        @type format: C{str}        @type format: C{str}
1254        @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.
1255        @rtype: C{str}        @rtype: C{str}
1256        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1257        """        """
1258        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1259           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1232  def whereNonZero(arg,tol=0.): Line 1287  def whereNonZero(arg,tol=0.):
1287    
1288     @param arg: argument     @param arg: argument
1289     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1290     @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.
1291     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1292     """     """
1293     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1294        if arg.rank==0:        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1295          if abs(arg)>tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1296             return numarray.array(1.)        return out
         else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
1297     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1298        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1299     elif isinstance(arg,float):     elif isinstance(arg,float):
1300        if abs(arg)>tol:        if abs(arg)>tol:
1301          return 1.          return 1.
# Line 1263  def whereNonZero(arg,tol=0.): Line 1311  def whereNonZero(arg,tol=0.):
1311     else:     else:
1312        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1313    
1314    def erf(arg):
1315       """
1316       returns erf of argument arg
1317    
1318       @param arg: argument
1319       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1320       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1321       @raises TypeError: if the type of the argument is not expected.
1322       """
1323       if isinstance(arg,escript.Data):
1324          return arg._erf()
1325       else:
1326          raise TypeError,"erf: Unknown argument type."
1327    
1328  def sin(arg):  def sin(arg):
1329     """     """
1330     returns sine of argument arg     returns sine of argument arg
1331    
1332     @param arg: argument     @param arg: argument
1333     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1334     @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.
1335     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1336     """     """
1337     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1307  class Sin_Symbol(DependendSymbol): Line 1369  class Sin_Symbol(DependendSymbol):
1369        @type format: C{str}        @type format: C{str}
1370        @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.
1371        @rtype: C{str}        @rtype: C{str}
1372        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1373        """        """
1374        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1375            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1359  def cos(arg): Line 1421  def cos(arg):
1421    
1422     @param arg: argument     @param arg: argument
1423     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1424     @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.
1425     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1426     """     """
1427     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1397  class Cos_Symbol(DependendSymbol): Line 1459  class Cos_Symbol(DependendSymbol):
1459        @type format: C{str}        @type format: C{str}
1460        @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.
1461        @rtype: C{str}        @rtype: C{str}
1462        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1463        """        """
1464        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1465            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1449  def tan(arg): Line 1511  def tan(arg):
1511    
1512     @param arg: argument     @param arg: argument
1513     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1514     @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.
1515     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1516     """     """
1517     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1487  class Tan_Symbol(DependendSymbol): Line 1549  class Tan_Symbol(DependendSymbol):
1549        @type format: C{str}        @type format: C{str}
1550        @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.
1551        @rtype: C{str}        @rtype: C{str}
1552        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1553        """        """
1554        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1555            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1539  def asin(arg): Line 1601  def asin(arg):
1601    
1602     @param arg: argument     @param arg: argument
1603     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1604     @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.
1605     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1606     """     """
1607     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1577  class Asin_Symbol(DependendSymbol): Line 1639  class Asin_Symbol(DependendSymbol):
1639        @type format: C{str}        @type format: C{str}
1640        @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.
1641        @rtype: C{str}        @rtype: C{str}
1642        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1643        """        """
1644        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1645            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1629  def acos(arg): Line 1691  def acos(arg):
1691    
1692     @param arg: argument     @param arg: argument
1693     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1694     @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.
1695     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1696     """     """
1697     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1667  class Acos_Symbol(DependendSymbol): Line 1729  class Acos_Symbol(DependendSymbol):
1729        @type format: C{str}        @type format: C{str}
1730        @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.
1731        @rtype: C{str}        @rtype: C{str}
1732        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1733        """        """
1734        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1735            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1719  def atan(arg): Line 1781  def atan(arg):
1781    
1782     @param arg: argument     @param arg: argument
1783     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1784     @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.
1785     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1786     """     """
1787     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1757  class Atan_Symbol(DependendSymbol): Line 1819  class Atan_Symbol(DependendSymbol):
1819        @type format: C{str}        @type format: C{str}
1820        @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.
1821        @rtype: C{str}        @rtype: C{str}
1822        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1823        """        """
1824        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1825            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1809  def sinh(arg): Line 1871  def sinh(arg):
1871    
1872     @param arg: argument     @param arg: argument
1873     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1874     @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.
1875     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1876     """     """
1877     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1847  class Sinh_Symbol(DependendSymbol): Line 1909  class Sinh_Symbol(DependendSymbol):
1909        @type format: C{str}        @type format: C{str}
1910        @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.
1911        @rtype: C{str}        @rtype: C{str}
1912        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1913        """        """
1914        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1915            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1899  def cosh(arg): Line 1961  def cosh(arg):
1961    
1962     @param arg: argument     @param arg: argument
1963     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1964     @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.
1965     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1966     """     """
1967     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1937  class Cosh_Symbol(DependendSymbol): Line 1999  class Cosh_Symbol(DependendSymbol):
1999        @type format: C{str}        @type format: C{str}
2000        @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.
2001        @rtype: C{str}        @rtype: C{str}
2002        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2003        """        """
2004        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2005            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1989  def tanh(arg): Line 2051  def tanh(arg):
2051    
2052     @param arg: argument     @param arg: argument
2053     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2054     @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.
2055     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2056     """     """
2057     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2027  class Tanh_Symbol(DependendSymbol): Line 2089  class Tanh_Symbol(DependendSymbol):
2089        @type format: C{str}        @type format: C{str}
2090        @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.
2091        @rtype: C{str}        @rtype: C{str}
2092        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2093        """        """
2094        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2095            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2079  def asinh(arg): Line 2141  def asinh(arg):
2141    
2142     @param arg: argument     @param arg: argument
2143     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2144     @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.
2145     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2146     """     """
2147     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2117  class Asinh_Symbol(DependendSymbol): Line 2179  class Asinh_Symbol(DependendSymbol):
2179        @type format: C{str}        @type format: C{str}
2180        @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.
2181        @rtype: C{str}        @rtype: C{str}
2182        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2183        """        """
2184        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2185            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2169  def acosh(arg): Line 2231  def acosh(arg):
2231    
2232     @param arg: argument     @param arg: argument
2233     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2234     @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.
2235     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2236     """     """
2237     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2207  class Acosh_Symbol(DependendSymbol): Line 2269  class Acosh_Symbol(DependendSymbol):
2269        @type format: C{str}        @type format: C{str}
2270        @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.
2271        @rtype: C{str}        @rtype: C{str}
2272        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2273        """        """
2274        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2275            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2259  def atanh(arg): Line 2321  def atanh(arg):
2321    
2322     @param arg: argument     @param arg: argument
2323     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2324     @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.
2325     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2326     """     """
2327     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2297  class Atanh_Symbol(DependendSymbol): Line 2359  class Atanh_Symbol(DependendSymbol):
2359        @type format: C{str}        @type format: C{str}
2360        @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.
2361        @rtype: C{str}        @rtype: C{str}
2362        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2363        """        """
2364        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2365            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2349  def exp(arg): Line 2411  def exp(arg):
2411    
2412     @param arg: argument     @param arg: argument
2413     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2414     @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.
2415     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2416     """     """
2417     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2387  class Exp_Symbol(DependendSymbol): Line 2449  class Exp_Symbol(DependendSymbol):
2449        @type format: C{str}        @type format: C{str}
2450        @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.
2451        @rtype: C{str}        @rtype: C{str}
2452        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2453        """        """
2454        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2455            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2439  def sqrt(arg): Line 2501  def sqrt(arg):
2501    
2502     @param arg: argument     @param arg: argument
2503     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2504     @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.
2505     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2506     """     """
2507     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2477  class Sqrt_Symbol(DependendSymbol): Line 2539  class Sqrt_Symbol(DependendSymbol):
2539        @type format: C{str}        @type format: C{str}
2540        @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.
2541        @rtype: C{str}        @rtype: C{str}
2542        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2543        """        """
2544        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2545            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2529  def log(arg): Line 2591  def log(arg):
2591    
2592     @param arg: argument     @param arg: argument
2593     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2594     @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.
2595     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2596     """     """
2597     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2567  class Log_Symbol(DependendSymbol): Line 2629  class Log_Symbol(DependendSymbol):
2629        @type format: C{str}        @type format: C{str}
2630        @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.
2631        @rtype: C{str}        @rtype: C{str}
2632        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2633        """        """
2634        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2635            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2619  def sign(arg): Line 2681  def sign(arg):
2681    
2682     @param arg: argument     @param arg: argument
2683     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2684     @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.
2685     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2686     """     """
2687     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2667  class Abs_Symbol(DependendSymbol): Line 2729  class Abs_Symbol(DependendSymbol):
2729        @type format: C{str}        @type format: C{str}
2730        @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.
2731        @rtype: C{str}        @rtype: C{str}
2732        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2733        """        """
2734        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2735            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2719  def minval(arg): Line 2781  def minval(arg):
2781    
2782     @param arg: argument     @param arg: argument
2783     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2784     @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.
2785     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2786     """     """
2787     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2760  class Minval_Symbol(DependendSymbol): Line 2822  class Minval_Symbol(DependendSymbol):
2822        @type format: C{str}        @type format: C{str}
2823        @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.
2824        @rtype: C{str}        @rtype: C{str}
2825        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2826        """        """
2827        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2828            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2796  def maxval(arg): Line 2858  def maxval(arg):
2858    
2859     @param arg: argument     @param arg: argument
2860     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2861     @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.
2862     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2863     """     """
2864     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2837  class Maxval_Symbol(DependendSymbol): Line 2899  class Maxval_Symbol(DependendSymbol):
2899        @type format: C{str}        @type format: C{str}
2900        @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.
2901        @rtype: C{str}        @rtype: C{str}
2902        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2903        """        """
2904        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2905            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2873  def length(arg): Line 2935  def length(arg):
2935    
2936     @param arg: argument     @param arg: argument
2937     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2938     @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.
2939     """     """
2940     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2941    
2942    def trace(arg,axis_offset=0):
2943       """
2944       returns the trace of arg which the sum of arg[k,k] over k.
2945    
2946       @param arg: argument
2947       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2948       @param axis_offset: C{axis_offset} to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
2949                      C{axis_offset} and axis_offset+1 must be equal.
2950       @type axis_offset: C{int}
2951       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2952       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2953       """
2954       if isinstance(arg,numarray.NumArray):
2955          sh=arg.shape
2956          if len(sh)<2:
2957            raise ValueError,"rank of argument must be greater than 1"
2958          if axis_offset<0 or axis_offset>len(sh)-2:
2959            raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
2960          s1=1
2961          for i in range(axis_offset): s1*=sh[i]
2962          s2=1
2963          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2964          if not sh[axis_offset] == sh[axis_offset+1]:
2965            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2966          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
2967          out=numarray.zeros([s1,s2],numarray.Float64)
2968          for i1 in range(s1):
2969            for i2 in range(s2):
2970                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
2971          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2972          return out
2973       elif isinstance(arg,escript.Data):
2974          if arg.getRank()<2:
2975            raise ValueError,"rank of argument must be greater than 1"
2976          if axis_offset<0 or axis_offset>arg.getRank()-2:
2977            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
2978          s=list(arg.getShape())        
2979          if not s[axis_offset] == s[axis_offset+1]:
2980            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2981          return arg._trace(axis_offset)
2982       elif isinstance(arg,float):
2983          raise TypeError,"illegal argument type float."
2984       elif isinstance(arg,int):
2985          raise TypeError,"illegal argument type int."
2986       elif isinstance(arg,Symbol):
2987          return Trace_Symbol(arg,axis_offset)
2988       else:
2989          raise TypeError,"Unknown argument type."
2990    
2991    class Trace_Symbol(DependendSymbol):
2992       """
2993       L{Symbol} representing the result of the trace function
2994       """
2995       def __init__(self,arg,axis_offset=0):
2996          """
2997          initialization of trace L{Symbol} with argument arg
2998          @param arg: argument of function
2999          @type arg: L{Symbol}.
3000          @param axis_offset: C{axis_offset} to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
3001                      C{axis_offset} and axis_offset+1 must be equal.
3002          @type axis_offset: C{int}
3003          """
3004          if arg.getRank()<2:
3005            raise ValueError,"rank of argument must be greater than 1"
3006          if axis_offset<0 or axis_offset>arg.getRank()-2:
3007            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3008          s=list(arg.getShape())        
3009          if not s[axis_offset] == s[axis_offset+1]:
3010            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3011          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3012    
3013       def getMyCode(self,argstrs,format="escript"):
3014          """
3015          returns a program code that can be used to evaluate the symbol.
3016    
3017          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3018          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3019          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3020          @type format: C{str}
3021          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3022          @rtype: C{str}
3023          @raise NotImplementedError: if the requested format is not available
3024          """
3025          if format=="escript" or format=="str"  or format=="text":
3026             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3027          else:
3028             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3029    
3030       def substitute(self,argvals):
3031          """
3032          assigns new values to symbols in the definition of the symbol.
3033          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3034    
3035          @param argvals: new values assigned to symbols
3036          @type argvals: C{dict} with keywords of type L{Symbol}.
3037          @return: result of the substitution process. Operations are executed as much as possible.
3038          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3039          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3040          """
3041          if argvals.has_key(self):
3042             arg=argvals[self]
3043             if self.isAppropriateValue(arg):
3044                return arg
3045             else:
3046                raise TypeError,"%s: new value is not appropriate."%str(self)
3047          else:
3048             arg=self.getSubstitutedArguments(argvals)
3049             return trace(arg[0],axis_offset=arg[1])
3050    
3051       def diff(self,arg):
3052          """
3053          differential of this object
3054    
3055          @param arg: the derivative is calculated with respect to arg
3056          @type arg: L{escript.Symbol}
3057          @return: derivative with respect to C{arg}
3058          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3059          """
3060          if arg==self:
3061             return identity(self.getShape())
3062          else:
3063             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3064    
3065    def transpose(arg,axis_offset=None):
3066       """
3067       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3068    
3069       @param arg: argument
3070       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3071       @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3072                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3073       @type axis_offset: C{int}
3074       @return: transpose of arg
3075       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3076       """
3077       if isinstance(arg,numarray.NumArray):
3078          if axis_offset==None: axis_offset=int(arg.rank/2)
3079          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3080       elif isinstance(arg,escript.Data):
3081          r=arg.getRank()
3082          if axis_offset==None: axis_offset=int(r/2)
3083          if axis_offset<0 or axis_offset>r:
3084            raise ValueError,"axis_offset must be between 0 and %s"%r
3085          return arg._transpose(axis_offset)
3086       elif isinstance(arg,float):
3087          if not ( axis_offset==0 or axis_offset==None):
3088            raise ValueError,"axis_offset must be 0 for float argument"
3089          return arg
3090       elif isinstance(arg,int):
3091          if not ( axis_offset==0 or axis_offset==None):
3092            raise ValueError,"axis_offset must be 0 for int argument"
3093          return float(arg)
3094       elif isinstance(arg,Symbol):
3095          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3096          return Transpose_Symbol(arg,axis_offset)
3097       else:
3098          raise TypeError,"Unknown argument type."
3099    
3100    class Transpose_Symbol(DependendSymbol):
3101       """
3102       L{Symbol} representing the result of the transpose function
3103       """
3104       def __init__(self,arg,axis_offset=None):
3105          """
3106          initialization of transpose L{Symbol} with argument arg
3107    
3108          @param arg: argument of function
3109          @type arg: L{Symbol}.
3110          @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3111                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3112          @type axis_offset: C{int}
3113          """
3114          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3115          if axis_offset<0 or axis_offset>arg.getRank():
3116            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3117          s=arg.getShape()
3118          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3119    
3120       def getMyCode(self,argstrs,format="escript"):
3121          """
3122          returns a program code that can be used to evaluate the symbol.
3123    
3124          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3125          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3126          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3127          @type format: C{str}
3128          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3129          @rtype: C{str}
3130          @raise NotImplementedError: if the requested format is not available
3131          """
3132          if format=="escript" or format=="str"  or format=="text":
3133             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3134          else:
3135             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3136    
3137       def substitute(self,argvals):
3138          """
3139          assigns new values to symbols in the definition of the symbol.
3140          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3141    
3142          @param argvals: new values assigned to symbols
3143          @type argvals: C{dict} with keywords of type L{Symbol}.
3144          @return: result of the substitution process. Operations are executed as much as possible.
3145          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3146          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3147          """
3148          if argvals.has_key(self):
3149             arg=argvals[self]
3150             if self.isAppropriateValue(arg):
3151                return arg
3152             else:
3153                raise TypeError,"%s: new value is not appropriate."%str(self)
3154          else:
3155             arg=self.getSubstitutedArguments(argvals)
3156             return transpose(arg[0],axis_offset=arg[1])
3157    
3158       def diff(self,arg):
3159          """
3160          differential of this object
3161    
3162          @param arg: the derivative is calculated with respect to arg
3163          @type arg: L{escript.Symbol}
3164          @return: derivative with respect to C{arg}
3165          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3166          """
3167          if arg==self:
3168             return identity(self.getShape())
3169          else:
3170             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3171    
3172    def swap_axes(arg,axis0=0,axis1=1):
3173       """
3174       returns the swap of arg by swaping the components axis0 and axis1
3175    
3176       @param arg: argument
3177       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3178       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3179       @type axis0: C{int}
3180       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3181       @type axis1: C{int}
3182       @return: C{arg} with swaped components
3183       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3184       """
3185       if axis0 > axis1:
3186          axis0,axis1=axis1,axis0
3187       if isinstance(arg,numarray.NumArray):
3188          return numarray.swapaxes(arg,axis0,axis1)
3189       elif isinstance(arg,escript.Data):
3190          return arg._swap_axes(axis0,axis1)
3191       elif isinstance(arg,float):
3192          raise TyepError,"float argument is not supported."
3193       elif isinstance(arg,int):
3194          raise TyepError,"int argument is not supported."
3195       elif isinstance(arg,Symbol):
3196          return SwapAxes_Symbol(arg,axis0,axis1)
3197       else:
3198          raise TypeError,"Unknown argument type."
3199    
3200    class SwapAxes_Symbol(DependendSymbol):
3201       """
3202       L{Symbol} representing the result of the swap function
3203       """
3204       def __init__(self,arg,axis0=0,axis1=1):
3205          """
3206          initialization of swap L{Symbol} with argument arg
3207    
3208          @param arg: argument
3209          @type arg: L{Symbol}.
3210          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3211          @type axis0: C{int}
3212          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3213          @type axis1: C{int}
3214          """
3215          if arg.getRank()<2:
3216             raise ValueError,"argument must have at least rank 2."
3217          if axis0<0 or axis0>arg.getRank()-1:
3218             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3219          if axis1<0 or axis1>arg.getRank()-1:
3220             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3221          if axis0 == axis1:
3222             raise ValueError,"axis indices must be different."
3223          if axis0 > axis1:
3224             axis0,axis1=axis1,axis0
3225          s=arg.getShape()
3226          s_out=[]
3227          for i in range(len(s)):
3228             if i == axis0:
3229                s_out.append(s[axis1])
3230             elif i == axis1:
3231                s_out.append(s[axis0])
3232             else:
3233                s_out.append(s[i])
3234          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3235    
3236       def getMyCode(self,argstrs,format="escript"):
3237          """
3238          returns a program code that can be used to evaluate the symbol.
3239    
3240          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3241          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3242          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3243          @type format: C{str}
3244          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3245          @rtype: C{str}
3246          @raise NotImplementedError: if the requested format is not available
3247          """
3248          if format=="escript" or format=="str"  or format=="text":
3249             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3250          else:
3251             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3252    
3253       def substitute(self,argvals):
3254          """
3255          assigns new values to symbols in the definition of the symbol.
3256          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3257    
3258          @param argvals: new values assigned to symbols
3259          @type argvals: C{dict} with keywords of type L{Symbol}.
3260          @return: result of the substitution process. Operations are executed as much as possible.
3261          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3262          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3263          """
3264          if argvals.has_key(self):
3265             arg=argvals[self]
3266             if self.isAppropriateValue(arg):
3267                return arg
3268             else:
3269                raise TypeError,"%s: new value is not appropriate."%str(self)
3270          else:
3271             arg=self.getSubstitutedArguments(argvals)
3272             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3273    
3274       def diff(self,arg):
3275          """
3276          differential of this object
3277    
3278          @param arg: the derivative is calculated with respect to arg
3279          @type arg: L{escript.Symbol}
3280          @return: derivative with respect to C{arg}
3281          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3282          """
3283          if arg==self:
3284             return identity(self.getShape())
3285          else:
3286             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3287    
3288    def symmetric(arg):
3289        """
3290        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3291    
3292        @param arg: square matrix. Must have rank 2 or 4 and be square.
3293        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3294        @return: symmetric part of arg
3295        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3296        """
3297        if isinstance(arg,numarray.NumArray):
3298          if arg.rank==2:
3299            if not (arg.shape[0]==arg.shape[1]):
3300               raise ValueError,"argument must be square."
3301          elif arg.rank==4:
3302            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3303               raise ValueError,"argument must be square."
3304          else:
3305            raise ValueError,"rank 2 or 4 is required."
3306          return (arg+transpose(arg))/2
3307        elif isinstance(arg,escript.Data):
3308          if arg.getRank()==2:
3309            if not (arg.getShape()[0]==arg.getShape()[1]):
3310               raise ValueError,"argument must be square."
3311            return arg._symmetric()
3312          elif arg.getRank()==4:
3313            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3314               raise ValueError,"argument must be square."
3315            return arg._symmetric()
3316          else:
3317            raise ValueError,"rank 2 or 4 is required."
3318        elif isinstance(arg,float):
3319          return arg
3320        elif isinstance(arg,int):
3321          return float(arg)
3322        elif isinstance(arg,Symbol):
3323          if arg.getRank()==2:
3324            if not (arg.getShape()[0]==arg.getShape()[1]):
3325               raise ValueError,"argument must be square."
3326          elif arg.getRank()==4:
3327            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3328               raise ValueError,"argument must be square."
3329          else:
3330            raise ValueError,"rank 2 or 4 is required."
3331          return (arg+transpose(arg))/2
3332        else:
3333          raise TypeError,"symmetric: Unknown argument type."
3334    
3335    def nonsymmetric(arg):
3336        """
3337        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3338    
3339        @param arg: square matrix. Must have rank 2 or 4 and be square.
3340        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3341        @return: nonsymmetric part of arg
3342        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3343        """
3344        if isinstance(arg,numarray.NumArray):
3345          if arg.rank==2:
3346            if not (arg.shape[0]==arg.shape[1]):
3347               raise ValueError,"nonsymmetric: argument must be square."
3348          elif arg.rank==4:
3349            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3350               raise ValueError,"nonsymmetric: argument must be square."
3351          else:
3352            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3353          return (arg-transpose(arg))/2
3354        elif isinstance(arg,escript.Data):
3355          if arg.getRank()==2:
3356            if not (arg.getShape()[0]==arg.getShape()[1]):
3357               raise ValueError,"argument must be square."
3358            return arg._nonsymmetric()
3359          elif arg.getRank()==4:
3360            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3361               raise ValueError,"argument must be square."
3362            return arg._nonsymmetric()
3363          else:
3364            raise ValueError,"rank 2 or 4 is required."
3365        elif isinstance(arg,float):
3366          return arg
3367        elif isinstance(arg,int):
3368          return float(arg)
3369        elif isinstance(arg,Symbol):
3370          if arg.getRank()==2:
3371            if not (arg.getShape()[0]==arg.getShape()[1]):
3372               raise ValueError,"nonsymmetric: argument must be square."
3373          elif arg.getRank()==4:
3374            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3375               raise ValueError,"nonsymmetric: argument must be square."
3376          else:
3377            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3378          return (arg-transpose(arg))/2
3379        else:
3380          raise TypeError,"nonsymmetric: Unknown argument type."
3381    
3382    def inverse(arg):
3383        """
3384        returns the inverse of the square matrix arg.
3385    
3386        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3387        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3388        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3389        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3390        @note: for L{escript.Data} objects the dimension is restricted to 3.
3391        """
3392        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3393        if isinstance(arg,numarray.NumArray):
3394          return numarray.linear_algebra.inverse(arg)
3395        elif isinstance(arg,escript.Data):
3396          return escript_inverse(arg)
3397        elif isinstance(arg,float):
3398          return 1./arg
3399        elif isinstance(arg,int):
3400          return 1./float(arg)
3401        elif isinstance(arg,Symbol):
3402          return Inverse_Symbol(arg)
3403        else:
3404          raise TypeError,"inverse: Unknown argument type."
3405    
3406    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3407          "arg is a Data objects!!!"
3408          if not arg.getRank()==2:
3409            raise ValueError,"escript_inverse: argument must have rank 2"
3410          s=arg.getShape()      
3411          if not s[0] == s[1]:
3412            raise ValueError,"escript_inverse: argument must be a square matrix."
3413          out=escript.Data(0.,s,arg.getFunctionSpace())
3414          if s[0]==1:
3415              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3416                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3417              out[0,0]=1./arg[0,0]
3418          elif s[0]==2:
3419              A11=arg[0,0]
3420              A12=arg[0,1]
3421              A21=arg[1,0]
3422              A22=arg[1,1]
3423              D = A11*A22-A12*A21
3424              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3425                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3426              D=1./D
3427              out[0,0]= A22*D
3428              out[1,0]=-A21*D
3429              out[0,1]=-A12*D
3430              out[1,1]= A11*D
3431          elif s[0]==3:
3432              A11=arg[0,0]
3433              A21=arg[1,0]
3434              A31=arg[2,0]
3435              A12=arg[0,1]
3436              A22=arg[1,1]
3437              A32=arg[2,1]
3438              A13=arg[0,2]
3439              A23=arg[1,2]
3440              A33=arg[2,2]
3441              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3442              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3443                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3444              D=1./D
3445              out[0,0]=(A22*A33-A23*A32)*D
3446              out[1,0]=(A31*A23-A21*A33)*D
3447              out[2,0]=(A21*A32-A31*A22)*D
3448              out[0,1]=(A13*A32-A12*A33)*D
3449              out[1,1]=(A11*A33-A31*A13)*D
3450              out[2,1]=(A12*A31-A11*A32)*D
3451              out[0,2]=(A12*A23-A13*A22)*D
3452              out[1,2]=(A13*A21-A11*A23)*D
3453              out[2,2]=(A11*A22-A12*A21)*D
3454          else:
3455             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3456          return out
3457    
3458    class Inverse_Symbol(DependendSymbol):
3459       """
3460       L{Symbol} representing the result of the inverse function
3461       """
3462       def __init__(self,arg):
3463          """
3464          initialization of inverse L{Symbol} with argument arg
3465          @param arg: argument of function
3466          @type arg: L{Symbol}.
3467          """
3468          if not arg.getRank()==2:
3469            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3470          s=arg.getShape()
3471          if not s[0] == s[1]:
3472            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3473          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3474    
3475       def getMyCode(self,argstrs,format="escript"):
3476          """
3477          returns a program code that can be used to evaluate the symbol.
3478    
3479          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3480          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3481          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3482          @type format: C{str}
3483          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3484          @rtype: C{str}
3485          @raise NotImplementedError: if the requested format is not available
3486          """
3487          if format=="escript" or format=="str"  or format=="text":
3488             return "inverse(%s)"%argstrs[0]
3489          else:
3490             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3491    
3492       def substitute(self,argvals):
3493          """
3494          assigns new values to symbols in the definition of the symbol.
3495          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3496    
3497          @param argvals: new values assigned to symbols
3498          @type argvals: C{dict} with keywords of type L{Symbol}.
3499          @return: result of the substitution process. Operations are executed as much as possible.
3500          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3501          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3502          """
3503          if argvals.has_key(self):
3504             arg=argvals[self]
3505             if self.isAppropriateValue(arg):
3506                return arg
3507             else:
3508                raise TypeError,"%s: new value is not appropriate."%str(self)
3509          else:
3510             arg=self.getSubstitutedArguments(argvals)
3511             return inverse(arg[0])
3512    
3513       def diff(self,arg):
3514          """
3515          differential of this object
3516    
3517          @param arg: the derivative is calculated with respect to arg
3518          @type arg: L{escript.Symbol}
3519          @return: derivative with respect to C{arg}
3520          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3521          """
3522          if arg==self:
3523             return identity(self.getShape())
3524          else:
3525             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3526    
3527    def eigenvalues(arg):
3528        """
3529        returns the eigenvalues of the square matrix arg.
3530    
3531        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3532                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3533        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3534        @return: the eigenvalues in increasing order.
3535        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3536        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3537        """
3538        if isinstance(arg,numarray.NumArray):
3539          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3540          out.sort()
3541          return out
3542        elif isinstance(arg,escript.Data):
3543          return arg._eigenvalues()
3544        elif isinstance(arg,Symbol):
3545          if not arg.getRank()==2:
3546            raise ValueError,"eigenvalues: argument must have rank 2"
3547          s=arg.getShape()      
3548          if not s[0] == s[1]:
3549            raise ValueError,"eigenvalues: argument must be a square matrix."
3550          if s[0]==1:
3551              return arg[0]
3552          elif s[0]==2:
3553              arg1=symmetric(arg)
3554              A11=arg1[0,0]
3555              A12=arg1[0,1]
3556              A22=arg1[1,1]
3557              trA=(A11+A22)/2.
3558              A11-=trA
3559              A22-=trA
3560              s=sqrt(A12**2-A11*A22)
3561              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3562          elif s[0]==3:
3563              arg1=symmetric(arg)
3564              A11=arg1[0,0]
3565              A12=arg1[0,1]
3566              A22=arg1[1,1]
3567              A13=arg1[0,2]
3568              A23=arg1[1,2]
3569              A33=arg1[2,2]
3570              trA=(A11+A22+A33)/3.
3571              A11-=trA
3572              A22-=trA
3573              A33-=trA
3574              A13_2=A13**2
3575              A23_2=A23**2
3576              A12_2=A12**2
3577              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3578              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3579              sq_p=sqrt(p/3.)
3580              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
3581              sq_p*=2.
3582              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3583               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3584               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3585              return trA+sq_p*f
3586          else:
3587             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3588        elif isinstance(arg,float):
3589          return arg
3590        elif isinstance(arg,int):
3591          return float(arg)
3592        else:
3593          raise TypeError,"eigenvalues: Unknown argument type."
3594    
3595    def eigenvalues_and_eigenvectors(arg):
3596        """
3597        returns the eigenvalues and eigenvectors of the square matrix arg.
3598    
3599        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3600                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3601        @type arg: L{escript.Data}
3602        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3603                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3604                 the eigenvector coresponding to the i-th eigenvalue.
3605        @rtype: L{tuple} of L{escript.Data}.
3606        @note: The dimension is restricted to 3.
3607        """
3608        if isinstance(arg,numarray.NumArray):
3609          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3610        elif isinstance(arg,escript.Data):
3611          return arg._eigenvalues_and_eigenvectors()
3612        elif isinstance(arg,Symbol):
3613          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3614        elif isinstance(arg,float):
3615          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3616        elif isinstance(arg,int):
3617          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3618        else:
3619          raise TypeError,"eigenvalues: Unknown argument type."
3620  #=======================================================  #=======================================================
3621  #  Binary operations:  #  Binary operations:
3622  #=======================================================  #=======================================================
# Line 2936  class Add_Symbol(DependendSymbol): Line 3676  class Add_Symbol(DependendSymbol):
3676        @type format: C{str}        @type format: C{str}
3677        @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.
3678        @rtype: C{str}        @rtype: C{str}
3679        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3680        """        """
3681        if format=="str" or format=="text":        if format=="str" or format=="text":
3682           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 2995  def mult(arg0,arg1): Line 3735  def mult(arg0,arg1):
3735         """         """
3736         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3737         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3738            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3739         else:         else:
3740            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3741                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3035  class Mult_Symbol(DependendSymbol): Line 3775  class Mult_Symbol(DependendSymbol):
3775        @type format: C{str}        @type format: C{str}
3776        @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.
3777        @rtype: C{str}        @rtype: C{str}
3778        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3779        """        """
3780        if format=="str" or format=="text":        if format=="str" or format=="text":
3781           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3095  def quotient(arg0,arg1): Line 3835  def quotient(arg0,arg1):
3835         """         """
3836         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3837         if testForZero(args[0]):         if testForZero(args[0]):
3838            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3839         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3840            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3841               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3140  class Quotient_Symbol(DependendSymbol): Line 3880  class Quotient_Symbol(DependendSymbol):
3880        @type format: C{str}        @type format: C{str}
3881        @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.
3882        @rtype: C{str}        @rtype: C{str}
3883        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3884        """        """
3885        if format=="str" or format=="text":        if format=="str" or format=="text":
3886           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3201  def power(arg0,arg1): Line 3941  def power(arg0,arg1):
3941         """         """
3942         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3943         if testForZero(args[0]):         if testForZero(args[0]):
3944            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3945         elif testForZero(args[1]):         elif testForZero(args[1]):
3946            return numarray.ones(args[0],numarray.Float)            return numarray.ones(pokeShape(args[1]),numarray.Float64)
3947         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3948            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3949         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):
# Line 3244  class Power_Symbol(DependendSymbol): Line 3984  class Power_Symbol(DependendSymbol):
3984        @type format: C{str}        @type format: C{str}
3985        @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.
3986        @rtype: C{str}        @rtype: C{str}
3987        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3988        """        """
3989        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
3990           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3304  def maximum(*args): Line 4044  def maximum(*args):
4044         if out==None:         if out==None:
4045            out=a            out=a
4046         else:         else:
4047            m=whereNegative(out-a)            diff=add(a,-out)
4048            out=m*a+(1.-m)*out            out=add(out,mult(wherePositive(diff),diff))
4049      return out      return out
4050        
4051  def minimum(*arg):  def minimum(*args):
4052      """      """
4053      the minimum over arguments args      the minimum over arguments args
4054    
# Line 3322  def minimum(*arg): Line 4062  def minimum(*arg):
4062         if out==None:         if out==None:
4063            out=a            out=a
4064         else:         else:
4065            m=whereNegative(out-a)            diff=add(a,-out)
4066            out=m*out+(1.-m)*a            out=add(out,mult(whereNegative(diff),diff))
4067      return out      return out
4068    
4069    def clip(arg,minval=None,maxval=None):
4070        """
4071        cuts the values of arg between minval and maxval
4072    
4073        @param arg: argument
4074        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4075        @param minval: lower range. If None no lower range is applied
4076        @type minval: C{float} or C{None}
4077        @param maxval: upper range. If None no upper range is applied
4078        @type maxval: C{float} or C{None}
4079        @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4080                 less then maxval are unchanged.
4081        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4082        @raise ValueError: if minval>maxval
4083        """
4084        if not minval==None and not maxval==None:
4085           if minval>maxval:
4086              raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4087        if minval == None:
4088            tmp=arg
4089        else:
4090            tmp=maximum(minval,arg)
4091        if maxval == None:
4092            return tmp
4093        else:
4094            return minimum(tmp,maxval)
4095    
4096        
4097  def inner(arg0,arg1):  def inner(arg0,arg1):
4098      """      """
# Line 3340  def inner(arg0,arg1): Line 4108  def inner(arg0,arg1):
4108      @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}
4109      @param arg1: second argument      @param arg1: second argument
4110      @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}
4111      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4112      @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
4113      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4114      """      """
# Line 3348  def inner(arg0,arg1): Line 4116  def inner(arg0,arg1):
4116      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4117      if not sh0==sh1:      if not sh0==sh1:
4118          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4119      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4120    
4121    def outer(arg0,arg1):
4122        """
4123        the outer product of the two argument:
4124    
4125        out[t,s]=arg0[t]*arg1[s]
4126    
4127        where
4128    
4129            - s runs through arg0.Shape
4130            - t runs through arg1.Shape
4131    
4132        @param arg0: first argument
4133        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4134        @param arg1: second argument
4135        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4136        @return: the outer product of arg0 and arg1 at each data point
4137        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4138        """
4139        return generalTensorProduct(arg0,arg1,axis_offset=0)
4140    
4141  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4142      """      """
4143        see L{matrix_mult}
4144        """
4145        return matrix_mult(arg0,arg1)
4146    
4147    def matrix_mult(arg0,arg1):
4148        """
4149      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4150    
4151      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4152    
4153            or      or
4154    
4155      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4156    
4157      The second dimension of arg0 and the length of arg1 must match.      The second dimension of arg0 and the first dimension of arg1 must match.
4158    
4159      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4160      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3376  def matrixmult(arg0,arg1): Line 4170  def matrixmult(arg0,arg1):
4170          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4171      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4172          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4173      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4174    
4175  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4176      """      """
4177      the outer product of the two argument:      use L{tensor_mult}
   
     out[t,s]=arg0[t]*arg1[s]  
   
     where s runs through arg0.Shape  
           t runs through arg1.Shape  
   
     @param arg0: first argument  
     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}  
     @param arg1: second argument  
     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}  
     @return: the outer product of arg0 and arg1 at each data point  
     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input  
4178      """      """
4179      return generalTensorProduct(arg0,arg1,offset=0)      return tensor_mult(arg0,arg1)
   
4180    
4181  def tensormult(arg0,arg1):  def tensor_mult(arg0,arg1):
4182      """      """
4183      the tensor product of the two argument:      the tensor product of the two argument:
   
4184            
4185      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4186    
4187      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4188    
4189                   or      or
4190    
4191      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4192    
# Line 3415  def tensormult(arg0,arg1): Line 4195  def tensormult(arg0,arg1):
4195    
4196      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]
4197                                
4198                   or      or
4199    
4200      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]
4201    
4202                   or      or
4203    
4204      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]
4205    
4206      In the first case the the second dimension of arg0 and the length of arg1 must match and        In the first case the the second dimension of arg0 and the last dimension of arg1 must match and  
4207      in the second case the two last dimensions of arg0 must match the shape of arg1.      in the second case the two last dimensions of arg0 must match the two first dimensions of arg1.
4208    
4209      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4210      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3436  def tensormult(arg0,arg1): Line 4216  def tensormult(arg0,arg1):
4216      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
4217      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4218      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4219         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4220      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4221         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4222      else:      else:
4223          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4224    
4225  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4226      """      """
4227      generalized tensor product      generalized tensor product
4228    
4229      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4230    
4231      where s runs through arg0.Shape[:arg0.Rank-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[offset:]  
4232    
4233      In the first case the the second dimension of arg0 and the length of arg1 must match and            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4234      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4235            - t runs through arg1.Shape[axis_offset:]
4236    
4237      @param arg0: first argument      @param arg0: first argument
4238      @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}
4239      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4240      @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}
4241      @return: the general tensor product of arg0 and arg1 at each data point.      @return: the general tensor product of arg0 and arg1 at each data point.
4242      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
# Line 3467  def generalTensorProduct(arg0,arg1,offse Line 4246  def generalTensorProduct(arg0,arg1,offse
4246      # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols      # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4247      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4248         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4249             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4250         else:         else:
4251             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4252                 raise ValueError,"generalTensorProduct: dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)                 raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4253             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4254             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4255             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4256             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4257             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4258             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4259             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4260             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4261             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4262             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4263             for i0 in range(d0):             for i0 in range(d0):
4264                      for i1 in range(d1):                      for i1 in range(d1):
4265                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4266             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4267             return out             return out
4268      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4269         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4270             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4271         else:         else:
4272             return escript_generalTensorProduct(arg0,arg1,offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,offset)             return escript_generalTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4273      else:            else:      
4274         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4275                                    
4276  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4277     """     """
4278     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4279     """     """
4280     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4281         """         """
4282         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4283    
4284         @param arg0: numerator         @param arg0: first argument
4285         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4286         @param arg1: denominator         @param arg1: second argument
4287         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4288         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4289         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4290         """         """
4291         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
4292         sh_arg1=pokeShape(arg1)         sh_arg1=pokeShape(arg1)
4293         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4294         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4295         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4296         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4297         if not sh01==sh10:         if not sh01==sh10:
4298             raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)             raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4299         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,offset])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4300    
4301     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4302        """        """
# Line 3529  class GeneralTensorProduct_Symbol(Depend Line 4308  class GeneralTensorProduct_Symbol(Depend
4308        @type format: C{str}        @type format: C{str}
4309        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4310        @rtype: C{str}        @rtype: C{str}
4311        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4312        """        """
4313        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4314           return "generalTensorProduct(%s,%s,offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4315        else:        else:
4316           raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)           raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4317    
# Line 3557  class GeneralTensorProduct_Symbol(Depend Line 4336  class GeneralTensorProduct_Symbol(Depend
4336           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4337           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4338    
4339  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4340      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4341      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4342      shape0=arg0.getShape()[:arg0.getRank()-offset]  
4343      shape01=arg0.getShape()[arg0.getRank()-offset:]  def transposed_matrix_mult(arg0,arg1):
4344      shape10=arg1.getShape()[:offset]      """
4345      shape1=arg1.getShape()[offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4346      if not shape01==shape10:  
4347          raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
4348    
4349      # create return value:      or
4350      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())  
4351      #      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4352      s0=[[]]  
4353      for k in shape0:      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4354            s=[]  
4355            for j in s0:      The first dimension of arg0 and arg1 must match.
4356                  for i in range(k): s.append(j+[slice(i,i)])  
4357            s0=s      @param arg0: first argument of rank 2
4358      s1=[[]]      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4359      for k in shape1:      @param arg1: second argument of at least rank 1
4360            s=[]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4361            for j in s1:      @return: the product of the transposed of arg0 and arg1 at each data point
4362                  for i in range(k): s.append(j+[slice(i,i)])      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4363            s1=s      @raise ValueError: if the shapes of the arguments are not appropriate
4364      s01=[[]]      """
4365      for k in shape01:      sh0=pokeShape(arg0)
4366            s=[]      sh1=pokeShape(arg1)
4367            for j in s01:      if not len(sh0)==2 :
4368                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"first argument must have rank 2"
4369            s01=s      if not len(sh1)==2 and not len(sh1)==1:
4370            raise ValueError,"second argument must have rank 1 or 2"
4371      for i0 in s0:      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4372         for i1 in s1:  
4373           s=escript.Scalar(0.,arg0.getFunctionSpace())  def transposed_tensor_mult(arg0,arg1):
4374           for i01 in s01:      """
4375              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      the tensor product of the transposed of the first and the second argument
4376           out.__setitem__(tuple(i0+i1),s)      
4377      return out      for arg0 of rank 2 this is
4378    
4379        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4380    
4381        or
4382    
4383        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4384    
4385      
4386        and for arg0 of rank 4 this is
4387    
4388        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4389                  
4390        or
4391    
4392        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4393    
4394        or
4395    
4396        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4397    
4398        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4399        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4400    
4401        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4402    
4403        @param arg0: first argument of rank 2 or 4
4404        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4405        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4406        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4407        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4408        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4409        """
4410        sh0=pokeShape(arg0)
4411        sh1=pokeShape(arg1)
4412        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4413           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4414        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4415           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4416        else:
4417            raise ValueError,"first argument must have rank 2 or 4"
4418    
4419    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4420        """
4421        generalized tensor product of transposed of arg0 and arg1:
4422    
4423        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4424    
4425        where
4426    
4427            - s runs through arg0.Shape[axis_offset:]
4428            - r runs trough arg0.Shape[:axis_offset]
4429            - t runs through arg1.Shape[axis_offset:]
4430    
4431        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4432        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4433    
4434        @param arg0: first argument
4435        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4436        @param arg1: second argument
4437        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4438        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4439        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4440        """
4441        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4442        arg0,arg1=matchType(arg0,arg1)
4443        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4444        if isinstance(arg0,numarray.NumArray):
4445           if isinstance(arg1,Symbol):
4446               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4447           else:
4448               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4449                   raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4450               arg0_c=arg0.copy()
4451               arg1_c=arg1.copy()
4452               sh0,sh1=arg0.shape,arg1.shape
4453               d0,d1,d01=1,1,1
4454               for i in sh0[axis_offset:]: d0*=i
4455               for i in sh1[axis_offset:]: d1*=i
4456               for i in sh0[:axis_offset]: d01*=i
4457               arg0_c.resize((d01,d0))
4458               arg1_c.resize((d01,d1))
4459               out=numarray.zeros((d0,d1),numarray.Float64)
4460               for i0 in range(d0):
4461                        for i1 in range(d1):
4462                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4463               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4464               return out
4465        elif isinstance(arg0,escript.Data):
4466           if isinstance(arg1,Symbol):
4467               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4468           else:
4469               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4470        else:      
4471           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4472                    
4473    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4474       """
4475       Symbol representing the general tensor product of the transposed of arg0 and arg1
4476       """
4477       def __init__(self,arg0,arg1,axis_offset=0):
4478           """
4479           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4480    
4481           @param arg0: first argument
4482           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4483           @param arg1: second argument
4484           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4485           @raise ValueError: inconsistent dimensions of arguments.
4486           @note: if both arguments have a spatial dimension, they must equal.
4487           """
4488           sh_arg0=pokeShape(arg0)
4489           sh_arg1=pokeShape(arg1)
4490           sh01=sh_arg0[:axis_offset]
4491           sh10=sh_arg1[:axis_offset]
4492           sh0=sh_arg0[axis_offset:]
4493           sh1=sh_arg1[axis_offset:]
4494           if not sh01==sh10:
4495               raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4496           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4497    
4498       def getMyCode(self,argstrs,format="escript"):
4499          """
4500          returns a program code that can be used to evaluate the symbol.
4501    
4502          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4503          @type argstrs: C{list} of length 2 of C{str}.
4504          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4505          @type format: C{str}
4506          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4507          @rtype: C{str}
4508          @raise NotImplementedError: if the requested format is not available
4509          """
4510          if format=="escript" or format=="str" or format=="text":
4511             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4512          else:
4513             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4514    
4515       def substitute(self,argvals):
4516          """
4517          assigns new values to symbols in the definition of the symbol.
4518          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4519    
4520          @param argvals: new values assigned to symbols
4521          @type argvals: C{dict} with keywords of type L{Symbol}.
4522          @return: result of the substitution process. Operations are executed as much as possible.
4523          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4524          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4525          """
4526          if argvals.has_key(self):
4527             arg=argvals[self]
4528             if self.isAppropriateValue(arg):
4529                return arg
4530             else:
4531                raise TypeError,"%s: new value is not appropriate."%str(self)
4532          else:
4533             args=self.getSubstitutedArguments(argvals)
4534             return generalTransposedTensorProduct(args[0],args[1],args[2])
4535    
4536    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4537        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4538        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4539    
4540    def matrix_transposed_mult(arg0,arg1):
4541        """
4542        matrix-transposed(matrix) product of the two argument:
4543    
4544        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4545    
4546        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4547    
4548        The last dimensions of arg0 and arg1 must match.
4549    
4550        @param arg0: first argument of rank 2
4551        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4552        @param arg1: second argument of rank 2
4553        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4554        @return: the product of arg0 and the transposed of arg1 at each data point
4555        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4556        @raise ValueError: if the shapes of the arguments are not appropriate
4557        """
4558        sh0=pokeShape(arg0)
4559        sh1=pokeShape(arg1)
4560        if not len(sh0)==2 :
4561            raise ValueError,"first argument must have rank 2"
4562        if not len(sh1)==2 and not len(sh1)==1:
4563            raise ValueError,"second argument must have rank 1 or 2"
4564        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4565    
4566    def tensor_transposed_mult(arg0,arg1):
4567        """
4568        the tensor product of the first and the transpose of the second argument
4569        
4570        for arg0 of rank 2 this is
4571    
4572        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4573    
4574        and for arg0 of rank 4 this is
4575    
4576        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4577                  
4578        or
4579    
4580        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4581    
4582        In the first case the the second dimension of arg0 and arg1 must match and  
4583        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4584    
4585        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4586    
4587        @param arg0: first argument of rank 2 or 4
4588        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4589        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4590        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4591        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4592        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4593        """
4594        sh0=pokeShape(arg0)
4595        sh1=pokeShape(arg1)
4596        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4597           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4598        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4599           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4600        else:
4601            raise ValueError,"first argument must have rank 2 or 4"
4602    
4603    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4604        """
4605        generalized tensor product of transposed of arg0 and arg1:
4606    
4607        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4608    
4609        where
4610    
4611            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4612            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4613            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4614    
4615        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4616        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4617    
4618        @param arg0: first argument
4619        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4620        @param arg1: second argument
4621        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4622        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4623        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4624        """
4625        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4626        arg0,arg1=matchType(arg0,arg1)
4627        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4628        if isinstance(arg0,numarray.NumArray):
4629           if isinstance(arg1,Symbol):
4630               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4631           else:
4632               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4633                   raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4634               arg0_c=arg0.copy()
4635               arg1_c=arg1.copy()
4636               sh0,sh1=arg0.shape,arg1.shape
4637               d0,d1,d01=1,1,1
4638               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4639               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4640               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4641               arg0_c.resize((d0,d01))
4642               arg1_c.resize((d1,d01))
4643               out=numarray.zeros((d0,d1),numarray.Float64)
4644               for i0 in range(d0):
4645                        for i1 in range(d1):
4646                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4647               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4648               return out
4649        elif isinstance(arg0,escript.Data):
4650           if isinstance(arg1,Symbol):
4651               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4652           else:
4653               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4654        else:      
4655           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4656                    
4657    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4658       """
4659       Symbol representing the general tensor product of arg0 and the transpose of arg1
4660       """
4661       def __init__(self,arg0,arg1,axis_offset=0):
4662           """
4663           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4664    
4665           @param arg0: first argument
4666           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4667           @param arg1: second argument
4668           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4669           @raise ValueError: inconsistent dimensions of arguments.
4670           @note: if both arguments have a spatial dimension, they must equal.
4671           """
4672           sh_arg0=pokeShape(arg0)
4673           sh_arg1=pokeShape(arg1)
4674           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4675           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4676           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4677           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4678           if not sh01==sh10:
4679               raise ValueError,"dimensions of last %s components in left argument don't match the last %s components in the right argument."%(axis_offset,axis_offset)
4680           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4681    
4682       def getMyCode(self,argstrs,format="escript"):
4683          """
4684          returns a program code that can be used to evaluate the symbol.
4685    
4686          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4687          @type argstrs: C{list} of length 2 of C{str}.
4688          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4689          @type format: C{str}
4690          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4691          @rtype: C{str}
4692          @raise NotImplementedError: if the requested format is not available
4693          """
4694          if format=="escript" or format=="str" or format=="text":
4695             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4696          else:
4697             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4698    
4699       def substitute(self,argvals):
4700          """
4701          assigns new values to symbols in the definition of the symbol.
4702          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4703    
4704          @param argvals: new values assigned to symbols
4705          @type argvals: C{dict} with keywords of type L{Symbol}.
4706          @return: result of the substitution process. Operations are executed as much as possible.
4707          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4708          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4709          """
4710          if argvals.has_key(self):
4711             arg=argvals[self]
4712             if self.isAppropriateValue(arg):
4713                return arg
4714             else:
4715                raise TypeError,"%s: new value is not appropriate."%str(self)
4716          else:
4717             args=self.getSubstitutedArguments(argvals)
4718             return generalTensorTransposedProduct(args[0],args[1],args[2])
4719    
4720    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4721        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4722        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4723    
4724  #=========================================================  #=========================================================
4725  #   some little helpers  #  functions dealing with spatial dependency
4726  #=========================================================  #=========================================================
4727  def grad(arg,where=None):  def grad(arg,where=None):
4728      """      """
4729      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4730    
4731        If C{g} is the returned object, then
4732    
4733      @param arg:   Data object representing the function which gradient        - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4734                    to be calculated.        - if C{arg} is rank 1 C{g[i,s]} is the derivative of C{arg[i]} with respect to the C{s}-th spatial dimension.
4735          - if C{arg} is rank 2 C{g[i,j,s]} is the derivative of C{arg[i,j]} with respect to the C{s}-th spatial dimension.
4736          - if C{arg} is rank 3 C{g[i,j,k,s]} is the derivative of C{arg[i,j,k]} with respect to the C{s}-th spatial dimension.
4737    
4738        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4739        @type arg: L{escript.Data} or L{Symbol}
4740      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4741                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4742        @type where: C{None} or L{escript.FunctionSpace}
4743        @return: gradient of arg.
4744        @rtype: L{escript.Data} or L{Symbol}
4745      """      """
4746      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4747         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3617  def grad(arg,where=None): Line 4751  def grad(arg,where=None):
4751         else:         else:
4752            return arg._grad(where)            return arg._grad(where)
4753      else:      else:
4754        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4755    
4756    class Grad_Symbol(DependendSymbol):
4757       """
4758       L{Symbol} representing the result of the gradient operator
4759       """
4760       def __init__(self,arg,where=None):
4761          """
4762          initialization of gradient L{Symbol} with argument arg
4763          @param arg: argument of function
4764          @type arg: L{Symbol}.
4765          @param where: FunctionSpace in which the gradient will be calculated.
4766                      If not present or C{None} an appropriate default is used.
4767          @type where: C{None} or L{escript.FunctionSpace}
4768          """
4769          d=arg.getDim()
4770          if d==None:
4771             raise ValueError,"argument must have a spatial dimension"
4772          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4773    
4774       def getMyCode(self,argstrs,format="escript"):
4775          """
4776          returns a program code that can be used to evaluate the symbol.
4777    
4778          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4779          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4780          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4781          @type format: C{str}
4782          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4783          @rtype: C{str}
4784          @raise NotImplementedError: if the requested format is not available
4785          """
4786          if format=="escript" or format=="str"  or format=="text":
4787             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4788          else:
4789             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4790    
4791       def substitute(self,argvals):
4792          """
4793          assigns new values to symbols in the definition of the symbol.
4794          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4795    
4796          @param argvals: new values assigned to symbols
4797          @type argvals: C{dict} with keywords of type L{Symbol}.
4798          @return: result of the substitution process. Operations are executed as much as possible.
4799          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4800          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4801          """
4802          if argvals.has_key(self):
4803             arg=argvals[self]
4804             if self.isAppropriateValue(arg):
4805                return arg
4806             else:
4807                raise TypeError,"%s: new value is not appropriate."%str(self)
4808          else:
4809             arg=self.getSubstitutedArguments(argvals)
4810             return grad(arg[0],where=arg[1])
4811    
4812       def diff(self,arg):
4813          """
4814          differential of this object
4815    
4816          @param arg: the derivative is calculated with respect to arg
4817          @type arg: L{escript.Symbol}
4818          @return: derivative with respect to C{arg}
4819          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4820          """
4821          if arg==self:
4822             return identity(self.getShape())
4823          else:
4824             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4825    
4826  def integrate(arg,where=None):  def integrate(arg,where=None):
4827      """      """
4828      Return the integral if the function represented by Data object arg over      Return the integral of the function C{arg} over its domain. If C{where} is present C{arg} is interpolated to C{where}
4829      its domain.      before integration.
4830    
4831      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4832        @type arg: L{escript.Data} or L{Symbol}
4833      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4834                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4835        @type where: C{None} or L{escript.FunctionSpace}
4836        @return: integral of arg.
4837        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4838      """      """
4839      if isinstance(arg,numarray.NumArray):      if isinstance(arg,Symbol):
         if checkForZero(arg):  
            return arg  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,float):  
         if checkForZero(arg):  
            return arg  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,int):  
         if checkForZero(arg):  
            return float(arg)  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,Symbol):  
4840         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4841      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4842         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
# Line 3654  def integrate(arg,where=None): Line 4847  def integrate(arg,where=None):
4847      else:      else:
4848        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4849    
4850    class Integrate_Symbol(DependendSymbol):
4851       """
4852       L{Symbol} representing the result of the spatial integration operator
4853       """
4854       def __init__(self,arg,where=None):
4855          """
4856          initialization of integration L{Symbol} with argument arg
4857          @param arg: argument of the integration
4858          @type arg: L{Symbol}.
4859          @param where: FunctionSpace in which the integration will be calculated.
4860                      If not present or C{None} an appropriate default is used.
4861          @type where: C{None} or L{escript.FunctionSpace}
4862          """
4863          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4864    
4865       def getMyCode(self,argstrs,format="escript"):
4866          """
4867          returns a program code that can be used to evaluate the symbol.
4868    
4869          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4870          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4871          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4872          @type format: C{str}
4873          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4874          @rtype: C{str}
4875          @raise NotImplementedError: if the requested format is not available
4876          """
4877          if format=="escript" or format=="str"  or format=="text":
4878             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4879          else:
4880             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4881    
4882       def substitute(self,argvals):
4883          """
4884          assigns new values to symbols in the definition of the symbol.
4885          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4886    
4887          @param argvals: new values assigned to symbols
4888          @type argvals: C{dict} with keywords of type L{Symbol}.
4889          @return: result of the substitution process. Operations are executed as much as possible.
4890          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4891          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4892          """
4893          if argvals.has_key(self):
4894             arg=argvals[self]
4895             if self.isAppropriateValue(arg):
4896                return arg
4897             else:
4898                raise TypeError,"%s: new value is not appropriate."%str(self)
4899          else:
4900             arg=self.getSubstitutedArguments(argvals)
4901             return integrate(arg[0],where=arg[1])
4902    
4903       def diff(self,arg):
4904          """
4905          differential of this object
4906    
4907          @param arg: the derivative is calculated with respect to arg
4908          @type arg: L{escript.Symbol}
4909          @return: derivative with respect to C{arg}
4910          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4911          """
4912          if arg==self:
4913             return identity(self.getShape())
4914          else:
4915             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4916    
4917    
4918  def interpolate(arg,where):  def interpolate(arg,where):
4919      """      """
4920      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
4921    
4922      @param arg:    interpolant      @param arg: interpolant
4923      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
4924        @param where: FunctionSpace to be interpolated to
4925        @type where: L{escript.FunctionSpace}
4926        @return: interpolated argument
4927        @rtype: C{escript.Data} or L{Symbol}
4928      """      """
4929      if testForZero(arg):      if isinstance(arg,Symbol):
4930        return 0         return Interpolate_Symbol(arg,where)
     elif isinstance(arg,Symbol):  
        return Interpolated_Symbol(arg,where)  
4931      else:      else:
4932         return escript.Data(arg,where)         return escript.Data(arg,where)
4933    
4934  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
4935      """     """
4936      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
4937       """
4938       def __init__(self,arg,where):
4939          """
4940          initialization of interpolation L{Symbol} with argument arg
4941          @param arg: argument of the interpolation
4942          @type arg: L{Symbol}.
4943          @param where: FunctionSpace into which the argument is interpolated.
4944          @type where: L{escript.FunctionSpace}
4945          """
4946          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4947    
4948      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
4949                    be calculated.        """
4950      @param where: FunctionSpace in which the gradient will be calculated.        returns a program code that can be used to evaluate the symbol.
                   If not present or C{None} an appropriate default is used.  
     """  
     g=grad(arg,where)  
     return trace(g,axis0=g.getRank()-2,axis1=g.getRank()-1)  
4951    
4952  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
4953      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4954      Returns the jump of arg across a continuity.        @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4955          @type format: C{str}
4956          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4957          @rtype: C{str}
4958          @raise NotImplementedError: if the requested format is not available
4959          """
4960          if format=="escript" or format=="str"  or format=="text":
4961             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
4962          else:
4963             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4964    
4965      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
4966                    to be calculated.        """
4967      """        assigns new values to symbols in the definition of the symbol.
4968      d=arg.getDomain()        The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
     return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())  
4969    
4970  #=============================        @param argvals: new values assigned to symbols
4971  #        @type argvals: C{dict} with keywords of type L{Symbol}.
4972  # wrapper for various functions: if the argument has attribute the function name        @return: result of the substitution process. Operations are executed as much as possible.
4973  # as an argument it calls the corresponding methods. Otherwise the corresponding        @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4974  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4975          """
4976          if argvals.has_key(self):
4977             arg=argvals[self]
4978             if self.isAppropriateValue(arg):
4979                return arg
4980             else:
4981                raise TypeError,"%s: new value is not appropriate."%str(self)
4982          else:
4983             arg=self.getSubstitutedArguments(argvals)
4984             return interpolate(arg[0],where=arg[1])
4985    
4986       def diff(self,arg):
4987          """
4988          differential of this object
4989    
4990          @param arg: the derivative is calculated with respect to arg
4991          @type arg: L{escript.Symbol}
4992          @return: derivative with respect to C{arg}
4993          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
4994          """
4995          if arg==self:
4996             return identity(self.getShape())
4997          else:
4998             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4999    
 # functions involving the underlying Domain:  
5000    
5001  def transpose(arg,axis=None):  def div(arg,where=None):
5002      """      """
5003      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
5004    
5005      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
5006        @type arg: L{escript.Data} or L{Symbol}
5007        @param where: FunctionSpace in which the divergence will be calculated.
5008                      If not present or C{None} an appropriate default is used.
5009        @type where: C{None} or L{escript.FunctionSpace}
5010        @return: divergence of arg.
5011        @rtype: L{escript.Data} or L{Symbol}
5012      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
5013      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5014         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
5015      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
5016         # hack for transpose          dim=arg.getDomain().getDim()
        r=arg.getRank()  
        if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"  
        s=arg.getShape()  
        out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())  
        for i in range(s[0]):  
           for j in range(s[1]):  
              out[j,i]=arg[i,j]  
        return out  
        # end hack for transpose  
        return arg.transpose(axis)  
5017      else:      else:
5018         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
5019        if not arg.getShape()==(dim,):
5020          raise ValueError,"div: expected shape is (%s,)"%dim
5021        return trace(grad(arg,where))
5022    
5023  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
5024      """      """
5025      Return      returns the jump of arg across the continuity of the domain
5026    
5027      @param arg:      @param arg: argument
5028        @type arg: L{escript.Data} or L{Symbol}
5029        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5030                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5031        @type domain: C{None} or L{escript.Domain}
5032        @return: jump of arg
5033        @rtype: L{escript.Data} or L{Symbol}
5034      """      """
5035      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
5036         s=list(arg.getShape())              return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
        s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])  
        return Trace_Symbol(arg,axis0=axis0,axis1=axis1)  
     elif isinstance(arg,escript.Data):  
        # hack for trace  
        s=arg.getShape()  
        if s[axis0]!=s[axis1]:  
            raise ValueError,"illegal axis in trace"  
        out=escript.Scalar(0.,arg.getFunctionSpace())  
        for i in range(s[axis0]):  
           out+=arg[i,i]  
        return out  
        # end hack for trace  
     else:  
        return numarray.trace(arg,axis0=axis0,axis1=axis1)  
5037    
5038    def L2(arg):
5039        """
5040        returns the L2 norm of arg at where
5041        
5042        @param arg: function which L2 to be calculated.
5043        @type arg: L{escript.Data} or L{Symbol}
5044        @return: L2 norm of arg.
5045        @rtype: L{float} or L{Symbol}
5046        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5047        """
5048        return sqrt(integrate(inner(arg,arg)))
5049    #=============================
5050    #
5051    
5052  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5053      """      """
5054      resorts the component of arg according to index      resorts the component of arg according to index
5055    
5056      """      """
5057      pass      raise NotImplementedError
5058  #  #
5059  # $Log: util.py,v $  # $Log: util.py,v $
5060  # Revision 1.14.2.16  2005/10/19 06:09:57  gross  # Revision 1.14.2.16  2005/10/19 06:09:57  gross

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

  ViewVC Help
Powered by ViewVC 1.1.26