/[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 876 by ksteube, Thu Oct 19 03:50:23 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=0.,maxval=1.):
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
4076        @type minval: C{float}
4077        @param maxval: upper range
4078        @type maxval: C{float}
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 minval>maxval:
4085           raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4086        return minimum(maximum(minval,arg),maxval)
4087    
4088        
4089  def inner(arg0,arg1):  def inner(arg0,arg1):
4090      """      """
# Line 3340  def inner(arg0,arg1): Line 4100  def inner(arg0,arg1):
4100      @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}
4101      @param arg1: second argument      @param arg1: second argument
4102      @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}
4103      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4104      @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
4105      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4106      """      """
# Line 3348  def inner(arg0,arg1): Line 4108  def inner(arg0,arg1):
4108      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4109      if not sh0==sh1:      if not sh0==sh1:
4110          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4111      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4112    
4113    def outer(arg0,arg1):
4114        """
4115        the outer product of the two argument:
4116    
4117        out[t,s]=arg0[t]*arg1[s]
4118    
4119        where
4120    
4121            - s runs through arg0.Shape
4122            - t runs through arg1.Shape
4123    
4124        @param arg0: first argument
4125        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4126        @param arg1: second argument
4127        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4128        @return: the outer product of arg0 and arg1 at each data point
4129        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4130        """
4131        return generalTensorProduct(arg0,arg1,axis_offset=0)
4132    
4133  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4134      """      """
4135        see L{matrix_mult}
4136        """
4137        return matrix_mult(arg0,arg1)
4138    
4139    def matrix_mult(arg0,arg1):
4140        """
4141      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4142    
4143      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4144    
4145            or      or
4146    
4147      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4148    
4149      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.
4150    
4151      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4152      @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 4162  def matrixmult(arg0,arg1):
4162          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4163      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4164          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4165      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4166    
4167  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4168      """      """
4169      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  
4170      """      """
4171      return generalTensorProduct(arg0,arg1,offset=0)      return tensor_mult(arg0,arg1)
4172    
4173    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4174      """      """
4175      the tensor product of the two argument:      the tensor product of the two argument:
   
4176            
4177      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4178    
4179      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4180    
4181                   or      or
4182    
4183      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4184    
# Line 3415  def tensormult(arg0,arg1): Line 4187  def tensormult(arg0,arg1):
4187    
4188      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]
4189                                
4190                   or      or
4191    
4192      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]
4193    
4194                   or      or
4195    
4196      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]
4197    
4198      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  
4199      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.
4200    
4201      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4202      @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 4208  def tensormult(arg0,arg1):
4208      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
4209      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4210      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4211         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4212      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):
4213         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4214      else:      else:
4215          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4216    
4217  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4218      """      """
4219      generalized tensor product      generalized tensor product
4220    
4221      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4222    
4223      where s runs through arg0.Shape[:arg0.Rank-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[offset:]  
4224    
4225      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]
4226      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4227            - t runs through arg1.Shape[axis_offset:]
4228    
4229      @param arg0: first argument      @param arg0: first argument
4230      @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}
4231      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4232      @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}
4233      @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.
4234      @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 4238  def generalTensorProduct(arg0,arg1,offse
4238      # 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
4239      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4240         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4241             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4242         else:         else:
4243             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4244                 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)
4245             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4246             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4247             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4248             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4249             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4250             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4251             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4252             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4253             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4254             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4255             for i0 in range(d0):             for i0 in range(d0):
4256                      for i1 in range(d1):                      for i1 in range(d1):
4257                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4258             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4259             return out             return out
4260      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4261         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4262             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4263         else:         else:
4264             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)
4265      else:            else:      
4266         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4267                                    
4268  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4269     """     """
4270     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4271     """     """
4272     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4273         """         """
4274         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4275    
4276         @param arg0: numerator         @param arg0: first argument
4277         @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}.
4278         @param arg1: denominator         @param arg1: second argument
4279         @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}.
4280         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4281         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4282         """         """
4283         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
4284         sh_arg1=pokeShape(arg1)         sh_arg1=pokeShape(arg1)
4285         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4286         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4287         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4288         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4289         if not sh01==sh10:         if not sh01==sh10:
4290             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)
4291         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])
4292    
4293     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4294        """        """
# Line 3529  class GeneralTensorProduct_Symbol(Depend Line 4300  class GeneralTensorProduct_Symbol(Depend
4300        @type format: C{str}        @type format: C{str}
4301        @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.
4302        @rtype: C{str}        @rtype: C{str}
4303        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4304        """        """
4305        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4306           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])
4307        else:        else:
4308           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)
4309    
# Line 3557  class GeneralTensorProduct_Symbol(Depend Line 4328  class GeneralTensorProduct_Symbol(Depend
4328           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4329           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4330    
4331  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4332      "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!!!"
4333      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4334      shape0=arg0.getShape()[:arg0.getRank()-offset]  
4335      shape01=arg0.getShape()[arg0.getRank()-offset:]  def transposed_matrix_mult(arg0,arg1):
4336      shape10=arg1.getShape()[:offset]      """
4337      shape1=arg1.getShape()[offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4338      if not shape01==shape10:  
4339          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]
4340    
4341      # create return value:      or
4342      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())  
4343      #      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4344      s0=[[]]  
4345      for k in shape0:      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4346            s=[]  
4347            for j in s0:      The first dimension of arg0 and arg1 must match.
4348                  for i in range(k): s.append(j+[slice(i,i)])  
4349            s0=s      @param arg0: first argument of rank 2
4350      s1=[[]]      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4351      for k in shape1:      @param arg1: second argument of at least rank 1
4352            s=[]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4353            for j in s1:      @return: the product of the transposed of arg0 and arg1 at each data point
4354                  for i in range(k): s.append(j+[slice(i,i)])      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4355            s1=s      @raise ValueError: if the shapes of the arguments are not appropriate
4356      s01=[[]]      """
4357      for k in shape01:      sh0=pokeShape(arg0)
4358            s=[]      sh1=pokeShape(arg1)
4359            for j in s01:      if not len(sh0)==2 :
4360                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"first argument must have rank 2"
4361            s01=s      if not len(sh1)==2 and not len(sh1)==1:
4362            raise ValueError,"second argument must have rank 1 or 2"
4363      for i0 in s0:      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4364         for i1 in s1:  
4365           s=escript.Scalar(0.,arg0.getFunctionSpace())  def transposed_tensor_mult(arg0,arg1):
4366           for i01 in s01:      """
4367              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      the tensor product of the transposed of the first and the second argument
4368           out.__setitem__(tuple(i0+i1),s)      
4369      return out      for arg0 of rank 2 this is
4370    
4371        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4372    
4373        or
4374    
4375        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4376    
4377      
4378        and for arg0 of rank 4 this is
4379    
4380        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4381                  
4382        or
4383    
4384        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4385    
4386        or
4387    
4388        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4389    
4390        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4391        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4392    
4393        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4394    
4395        @param arg0: first argument of rank 2 or 4
4396        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4397        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4398        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4399        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4400        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4401        """
4402        sh0=pokeShape(arg0)
4403        sh1=pokeShape(arg1)
4404        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4405           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4406        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4407           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4408        else:
4409            raise ValueError,"first argument must have rank 2 or 4"
4410    
4411    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4412        """
4413        generalized tensor product of transposed of arg0 and arg1:
4414    
4415        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4416    
4417        where
4418    
4419            - s runs through arg0.Shape[axis_offset:]
4420            - r runs trough arg0.Shape[:axis_offset]
4421            - t runs through arg1.Shape[axis_offset:]
4422    
4423        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4424        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4425    
4426        @param arg0: first argument
4427        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4428        @param arg1: second argument
4429        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4430        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4431        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4432        """
4433        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4434        arg0,arg1=matchType(arg0,arg1)
4435        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4436        if isinstance(arg0,numarray.NumArray):
4437           if isinstance(arg1,Symbol):
4438               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4439           else:
4440               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4441                   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)
4442               arg0_c=arg0.copy()
4443               arg1_c=arg1.copy()
4444               sh0,sh1=arg0.shape,arg1.shape
4445               d0,d1,d01=1,1,1
4446               for i in sh0[axis_offset:]: d0*=i
4447               for i in sh1[axis_offset:]: d1*=i
4448               for i in sh0[:axis_offset]: d01*=i
4449               arg0_c.resize((d01,d0))
4450               arg1_c.resize((d01,d1))
4451               out=numarray.zeros((d0,d1),numarray.Float64)
4452               for i0 in range(d0):
4453                        for i1 in range(d1):
4454                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4455               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4456               return out
4457        elif isinstance(arg0,escript.Data):
4458           if isinstance(arg1,Symbol):
4459               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4460           else:
4461               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4462        else:      
4463           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4464                    
4465    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4466       """
4467       Symbol representing the general tensor product of the transposed of arg0 and arg1
4468       """
4469       def __init__(self,arg0,arg1,axis_offset=0):
4470           """
4471           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4472    
4473           @param arg0: first argument
4474           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4475           @param arg1: second argument
4476           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4477           @raise ValueError: inconsistent dimensions of arguments.
4478           @note: if both arguments have a spatial dimension, they must equal.
4479           """
4480           sh_arg0=pokeShape(arg0)
4481           sh_arg1=pokeShape(arg1)
4482           sh01=sh_arg0[:axis_offset]
4483           sh10=sh_arg1[:axis_offset]
4484           sh0=sh_arg0[axis_offset:]
4485           sh1=sh_arg1[axis_offset:]
4486           if not sh01==sh10:
4487               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)
4488           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4489    
4490       def getMyCode(self,argstrs,format="escript"):
4491          """
4492          returns a program code that can be used to evaluate the symbol.
4493    
4494          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4495          @type argstrs: C{list} of length 2 of C{str}.
4496          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4497          @type format: C{str}
4498          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4499          @rtype: C{str}
4500          @raise NotImplementedError: if the requested format is not available
4501          """
4502          if format=="escript" or format=="str" or format=="text":
4503             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4504          else:
4505             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4506    
4507       def substitute(self,argvals):
4508          """
4509          assigns new values to symbols in the definition of the symbol.
4510          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4511    
4512          @param argvals: new values assigned to symbols
4513          @type argvals: C{dict} with keywords of type L{Symbol}.
4514          @return: result of the substitution process. Operations are executed as much as possible.
4515          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4516          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4517          """
4518          if argvals.has_key(self):
4519             arg=argvals[self]
4520             if self.isAppropriateValue(arg):
4521                return arg
4522             else:
4523                raise TypeError,"%s: new value is not appropriate."%str(self)
4524          else:
4525             args=self.getSubstitutedArguments(argvals)
4526             return generalTransposedTensorProduct(args[0],args[1],args[2])
4527    
4528    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4529        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4530        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4531    
4532    def matrix_transposed_mult(arg0,arg1):
4533        """
4534        matrix-transposed(matrix) product of the two argument:
4535    
4536        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4537    
4538        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4539    
4540        The last dimensions of arg0 and arg1 must match.
4541    
4542        @param arg0: first argument of rank 2
4543        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4544        @param arg1: second argument of rank 2
4545        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4546        @return: the product of arg0 and the transposed of arg1 at each data point
4547        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4548        @raise ValueError: if the shapes of the arguments are not appropriate
4549        """
4550        sh0=pokeShape(arg0)
4551        sh1=pokeShape(arg1)
4552        if not len(sh0)==2 :
4553            raise ValueError,"first argument must have rank 2"
4554        if not len(sh1)==2 and not len(sh1)==1:
4555            raise ValueError,"second argument must have rank 1 or 2"
4556        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4557    
4558    def tensor_transposed_mult(arg0,arg1):
4559        """
4560        the tensor product of the first and the transpose of the second argument
4561        
4562        for arg0 of rank 2 this is
4563    
4564        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4565    
4566        and for arg0 of rank 4 this is
4567    
4568        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4569                  
4570        or
4571    
4572        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4573    
4574        In the first case the the second dimension of arg0 and arg1 must match and  
4575        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4576    
4577        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4578    
4579        @param arg0: first argument of rank 2 or 4
4580        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4581        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4582        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4583        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4584        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4585        """
4586        sh0=pokeShape(arg0)
4587        sh1=pokeShape(arg1)
4588        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4589           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4590        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4591           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4592        else:
4593            raise ValueError,"first argument must have rank 2 or 4"
4594    
4595    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4596        """
4597        generalized tensor product of transposed of arg0 and arg1:
4598    
4599        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4600    
4601        where
4602    
4603            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4604            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4605            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4606    
4607        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4608        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4609    
4610        @param arg0: first argument
4611        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4612        @param arg1: second argument
4613        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4614        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4615        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4616        """
4617        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4618        arg0,arg1=matchType(arg0,arg1)
4619        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4620        if isinstance(arg0,numarray.NumArray):
4621           if isinstance(arg1,Symbol):
4622               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4623           else:
4624               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4625                   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)
4626               arg0_c=arg0.copy()
4627               arg1_c=arg1.copy()
4628               sh0,sh1=arg0.shape,arg1.shape
4629               d0,d1,d01=1,1,1
4630               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4631               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4632               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4633               arg0_c.resize((d0,d01))
4634               arg1_c.resize((d1,d01))
4635               out=numarray.zeros((d0,d1),numarray.Float64)
4636               for i0 in range(d0):
4637                        for i1 in range(d1):
4638                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4639               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4640               return out
4641        elif isinstance(arg0,escript.Data):
4642           if isinstance(arg1,Symbol):
4643               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4644           else:
4645               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4646        else:      
4647           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4648                    
4649    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4650       """
4651       Symbol representing the general tensor product of arg0 and the transpose of arg1
4652       """
4653       def __init__(self,arg0,arg1,axis_offset=0):
4654           """
4655           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4656    
4657           @param arg0: first argument
4658           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4659           @param arg1: second argument
4660           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4661           @raise ValueError: inconsistent dimensions of arguments.
4662           @note: if both arguments have a spatial dimension, they must equal.
4663           """
4664           sh_arg0=pokeShape(arg0)
4665           sh_arg1=pokeShape(arg1)
4666           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4667           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4668           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4669           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4670           if not sh01==sh10:
4671               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)
4672           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4673    
4674       def getMyCode(self,argstrs,format="escript"):
4675          """
4676          returns a program code that can be used to evaluate the symbol.
4677    
4678          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4679          @type argstrs: C{list} of length 2 of C{str}.
4680          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4681          @type format: C{str}
4682          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4683          @rtype: C{str}
4684          @raise NotImplementedError: if the requested format is not available
4685          """
4686          if format=="escript" or format=="str" or format=="text":
4687             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4688          else:
4689             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4690    
4691       def substitute(self,argvals):
4692          """
4693          assigns new values to symbols in the definition of the symbol.
4694          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4695    
4696          @param argvals: new values assigned to symbols
4697          @type argvals: C{dict} with keywords of type L{Symbol}.
4698          @return: result of the substitution process. Operations are executed as much as possible.
4699          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4700          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4701          """
4702          if argvals.has_key(self):
4703             arg=argvals[self]
4704             if self.isAppropriateValue(arg):
4705                return arg
4706             else:
4707                raise TypeError,"%s: new value is not appropriate."%str(self)
4708          else:
4709             args=self.getSubstitutedArguments(argvals)
4710             return generalTensorTransposedProduct(args[0],args[1],args[2])
4711    
4712    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4713        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4714        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4715    
4716  #=========================================================  #=========================================================
4717  #   some little helpers  #  functions dealing with spatial dependency
4718  #=========================================================  #=========================================================
4719  def grad(arg,where=None):  def grad(arg,where=None):
4720      """      """
4721      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4722    
4723        If C{g} is the returned object, then
4724    
4725      @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.
4726                    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.
4727          - 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.
4728          - 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.
4729    
4730        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4731        @type arg: L{escript.Data} or L{Symbol}
4732      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4733                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4734        @type where: C{None} or L{escript.FunctionSpace}
4735        @return: gradient of arg.
4736        @rtype: L{escript.Data} or L{Symbol}
4737      """      """
4738      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4739         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3617  def grad(arg,where=None): Line 4743  def grad(arg,where=None):
4743         else:         else:
4744            return arg._grad(where)            return arg._grad(where)
4745      else:      else:
4746        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4747    
4748    class Grad_Symbol(DependendSymbol):
4749       """
4750       L{Symbol} representing the result of the gradient operator
4751       """
4752       def __init__(self,arg,where=None):
4753          """
4754          initialization of gradient L{Symbol} with argument arg
4755          @param arg: argument of function
4756          @type arg: L{Symbol}.
4757          @param where: FunctionSpace in which the gradient will be calculated.
4758                      If not present or C{None} an appropriate default is used.
4759          @type where: C{None} or L{escript.FunctionSpace}
4760          """
4761          d=arg.getDim()
4762          if d==None:
4763             raise ValueError,"argument must have a spatial dimension"
4764          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4765    
4766       def getMyCode(self,argstrs,format="escript"):
4767          """
4768          returns a program code that can be used to evaluate the symbol.
4769    
4770          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4771          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4772          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4773          @type format: C{str}
4774          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4775          @rtype: C{str}
4776          @raise NotImplementedError: if the requested format is not available
4777          """
4778          if format=="escript" or format=="str"  or format=="text":
4779             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4780          else:
4781             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4782    
4783       def substitute(self,argvals):
4784          """
4785          assigns new values to symbols in the definition of the symbol.
4786          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4787    
4788          @param argvals: new values assigned to symbols
4789          @type argvals: C{dict} with keywords of type L{Symbol}.
4790          @return: result of the substitution process. Operations are executed as much as possible.
4791          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4792          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4793          """
4794          if argvals.has_key(self):
4795             arg=argvals[self]
4796             if self.isAppropriateValue(arg):
4797                return arg
4798             else:
4799                raise TypeError,"%s: new value is not appropriate."%str(self)
4800          else:
4801             arg=self.getSubstitutedArguments(argvals)
4802             return grad(arg[0],where=arg[1])
4803    
4804       def diff(self,arg):
4805          """
4806          differential of this object
4807    
4808          @param arg: the derivative is calculated with respect to arg
4809          @type arg: L{escript.Symbol}
4810          @return: derivative with respect to C{arg}
4811          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4812          """
4813          if arg==self:
4814             return identity(self.getShape())
4815          else:
4816             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4817    
4818  def integrate(arg,where=None):  def integrate(arg,where=None):
4819      """      """
4820      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}
4821      its domain.      before integration.
4822    
4823      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4824        @type arg: L{escript.Data} or L{Symbol}
4825      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4826                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4827        @type where: C{None} or L{escript.FunctionSpace}
4828        @return: integral of arg.
4829        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4830      """      """
4831      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):  
4832         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4833      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4834         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 4839  def integrate(arg,where=None):
4839      else:      else:
4840        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4841    
4842    class Integrate_Symbol(DependendSymbol):
4843       """
4844       L{Symbol} representing the result of the spatial integration operator
4845       """
4846       def __init__(self,arg,where=None):
4847          """
4848          initialization of integration L{Symbol} with argument arg
4849          @param arg: argument of the integration
4850          @type arg: L{Symbol}.
4851          @param where: FunctionSpace in which the integration will be calculated.
4852                      If not present or C{None} an appropriate default is used.
4853          @type where: C{None} or L{escript.FunctionSpace}
4854          """
4855          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4856    
4857       def getMyCode(self,argstrs,format="escript"):
4858          """
4859          returns a program code that can be used to evaluate the symbol.
4860    
4861          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4862          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4863          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4864          @type format: C{str}
4865          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4866          @rtype: C{str}
4867          @raise NotImplementedError: if the requested format is not available
4868          """
4869          if format=="escript" or format=="str"  or format=="text":
4870             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4871          else:
4872             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4873    
4874       def substitute(self,argvals):
4875          """
4876          assigns new values to symbols in the definition of the symbol.
4877          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4878    
4879          @param argvals: new values assigned to symbols
4880          @type argvals: C{dict} with keywords of type L{Symbol}.
4881          @return: result of the substitution process. Operations are executed as much as possible.
4882          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4883          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4884          """
4885          if argvals.has_key(self):
4886             arg=argvals[self]
4887             if self.isAppropriateValue(arg):
4888                return arg
4889             else:
4890                raise TypeError,"%s: new value is not appropriate."%str(self)
4891          else:
4892             arg=self.getSubstitutedArguments(argvals)
4893             return integrate(arg[0],where=arg[1])
4894    
4895       def diff(self,arg):
4896          """
4897          differential of this object
4898    
4899          @param arg: the derivative is calculated with respect to arg
4900          @type arg: L{escript.Symbol}
4901          @return: derivative with respect to C{arg}
4902          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4903          """
4904          if arg==self:
4905             return identity(self.getShape())
4906          else:
4907             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4908    
4909    
4910  def interpolate(arg,where):  def interpolate(arg,where):
4911      """      """
4912      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
4913    
4914      @param arg:    interpolant      @param arg: interpolant
4915      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
4916        @param where: FunctionSpace to be interpolated to
4917        @type where: L{escript.FunctionSpace}
4918        @return: interpolated argument
4919        @rtype: C{escript.Data} or L{Symbol}
4920      """      """
4921      if testForZero(arg):      if isinstance(arg,Symbol):
4922        return 0         return Interpolate_Symbol(arg,where)
     elif isinstance(arg,Symbol):  
        return Interpolated_Symbol(arg,where)  
4923      else:      else:
4924         return escript.Data(arg,where)         return escript.Data(arg,where)
4925    
4926  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
4927      """     """
4928      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
4929       """
4930       def __init__(self,arg,where):
4931          """
4932          initialization of interpolation L{Symbol} with argument arg
4933          @param arg: argument of the interpolation
4934          @type arg: L{Symbol}.
4935          @param where: FunctionSpace into which the argument is interpolated.
4936          @type where: L{escript.FunctionSpace}
4937          """
4938          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4939    
4940      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
4941                    be calculated.        """
4942      @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)  
4943    
4944  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
4945      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4946      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.
4947          @type format: C{str}
4948          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4949          @rtype: C{str}
4950          @raise NotImplementedError: if the requested format is not available
4951          """
4952          if format=="escript" or format=="str"  or format=="text":
4953             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
4954          else:
4955             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4956    
4957      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
4958                    to be calculated.        """
4959      """        assigns new values to symbols in the definition of the symbol.
4960      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())  
4961    
4962  #=============================        @param argvals: new values assigned to symbols
4963  #        @type argvals: C{dict} with keywords of type L{Symbol}.
4964  # 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.
4965  # 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
4966  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4967          """
4968          if argvals.has_key(self):
4969             arg=argvals[self]
4970             if self.isAppropriateValue(arg):
4971                return arg
4972             else:
4973                raise TypeError,"%s: new value is not appropriate."%str(self)
4974          else:
4975             arg=self.getSubstitutedArguments(argvals)
4976             return interpolate(arg[0],where=arg[1])
4977    
4978       def diff(self,arg):
4979          """
4980          differential of this object
4981    
4982          @param arg: the derivative is calculated with respect to arg
4983          @type arg: L{escript.Symbol}
4984          @return: derivative with respect to C{arg}
4985          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
4986          """
4987          if arg==self:
4988             return identity(self.getShape())
4989          else:
4990             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4991    
 # functions involving the underlying Domain:  
4992    
4993  def transpose(arg,axis=None):  def div(arg,where=None):
4994      """      """
4995      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
4996    
4997      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
4998        @type arg: L{escript.Data} or L{Symbol}
4999        @param where: FunctionSpace in which the divergence will be calculated.
5000                      If not present or C{None} an appropriate default is used.
5001        @type where: C{None} or L{escript.FunctionSpace}
5002        @return: divergence of arg.
5003        @rtype: L{escript.Data} or L{Symbol}
5004      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
5005      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5006         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
5007      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
5008         # 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)  
5009      else:      else:
5010         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
5011        if not arg.getShape()==(dim,):
5012          raise ValueError,"div: expected shape is (%s,)"%dim
5013        return trace(grad(arg,where))
5014    
5015  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
5016      """      """
5017      Return      returns the jump of arg across the continuity of the domain
5018    
5019      @param arg:      @param arg: argument
5020        @type arg: L{escript.Data} or L{Symbol}
5021        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5022                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5023        @type domain: C{None} or L{escript.Domain}
5024        @return: jump of arg
5025        @rtype: L{escript.Data} or L{Symbol}
5026      """      """
5027      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
5028         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)  
5029    
5030    def L2(arg):
5031        """
5032        returns the L2 norm of arg at where
5033        
5034        @param arg: function which L2 to be calculated.
5035        @type arg: L{escript.Data} or L{Symbol}
5036        @return: L2 norm of arg.
5037        @rtype: L{float} or L{Symbol}
5038        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5039        """
5040        return sqrt(integrate(inner(arg,arg)))
5041    #=============================
5042    #
5043    
5044  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5045      """      """
5046      resorts the component of arg according to index      resorts the component of arg according to index
5047    
5048      """      """
5049      pass      raise NotImplementedError
5050  #  #
5051  # $Log: util.py,v $  # $Log: util.py,v $
5052  # 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.876

  ViewVC Help
Powered by ViewVC 1.1.26