/[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 353 by gross, Wed Dec 14 04:13:59 2005 UTC revision 1042 by gross, Mon Mar 19 03:50:34 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2    
3  """  """
4  Utility functions for escript  Utility functions for escript
5    
 @remark:  This module is under construction and is still tested!!!  
   
6  @var __author__: name of author  @var __author__: name of author
7  @var __licence__: licence agreement  @var __copyright__: copyrights
8    @var __license__: licence agreement
9  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
10  @var __version__: version  @var __version__: version
11  @var __date__: date of the version  @var __date__: date of the version
12  """  """
13                                                                                                                                                                                                                                                                                                                                                                                                            
14  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
15  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                        http://www.access.edu.au
17                    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19                 http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
21  __version__="$Revision$"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
# Line 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:
33  #=========================================================  #=========================================================
34    def getTagNames(domain):
35        """
36        returns a list of the tag names used by the domain
37    
38        
39        @param domain: a domain object
40        @type domain: C{escript.Domain}
41        @return: a list of the tag name used by the domain.
42        @rtype: C{list} of C{str}
43        """
44        return [n.strip() for n in domain.showTagNames().split(",") ]
45        
46  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
47      """      """
48      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.
49    
50      Example:      Example::
51    
52         tmp=Scalar(..)         tmp=Scalar(..)
53         v=Vector(..)         v=Vector(..)
# Line 96  def saveDX(filename,domain=None,**data): Line 75  def saveDX(filename,domain=None,**data):
75      """      """
76      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.
77    
78      Example:      Example::
79    
80         tmp=Scalar(..)         tmp=Scalar(..)
81         v=Vector(..)         v=Vector(..)
# Line 125  def kronecker(d=3): Line 104  def kronecker(d=3):
104     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
105    
106     @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
107     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
108     @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
109     @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}  
110     """     """
111     return identityTensor(d)     return identityTensor(d)
112    
# Line 143  def identity(shape=()): Line 121  def identity(shape=()):
121     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
122     """     """
123     if len(shape)>0:     if len(shape)>0:
124        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
125        if len(shape)==1:        if len(shape)==1:
126            for i0 in range(shape[0]):            for i0 in range(shape[0]):
127               out[i0,i0]=1.               out[i0,i0]=1.
   
128        elif len(shape)==2:        elif len(shape)==2:
129            for i0 in range(shape[0]):            for i0 in range(shape[0]):
130               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 163  def identityTensor(d=3): Line 140  def identityTensor(d=3):
140     return the dxd identity matrix     return the dxd identity matrix
141    
142     @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
143     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
144     @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
145     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
146     """     """
147     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
148        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
149     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
150           return identity((d.getDim(),))
151       else:
152           return identity((d,))
153    
154  def identityTensor4(d=3):  def identityTensor4(d=3):
155     """     """
# Line 178  def identityTensor4(d=3): Line 158  def identityTensor4(d=3):
158     @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
159     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
160     @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
161     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
162     """     """
163     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
164        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
165     return identity((d,d))     elif isinstance(d,escript.Domain):
166           return identity((d.getDim(),d.getDim()))
167       else:
168           return identity((d,d))
169    
170  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
171     """     """
# Line 191  def unitVector(i=0,d=3): Line 174  def unitVector(i=0,d=3):
174     @param i: index     @param i: index
175     @type i: C{int}     @type i: C{int}
176     @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
177     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
178     @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
179     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
180     """     """
181     return kronecker(d)[i]     return kronecker(d)[i]
182    
# Line 249  def inf(arg): Line 232  def inf(arg):
232    
233      @param arg: argument      @param arg: argument
234      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
235      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
236      @rtype: C{float}      @rtype: C{float}
237      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
238      """      """
# Line 268  def inf(arg): Line 251  def inf(arg):
251  #=========================================================================  #=========================================================================
252  #   some little helpers  #   some little helpers
253  #=========================================================================  #=========================================================================
254  def pokeShape(arg):  def getRank(arg):
255        """
256        identifies the rank of its argument
257    
258        @param arg: a given object
259        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
260        @return: the rank of the argument
261        @rtype: C{int}
262        @raise TypeError: if type of arg cannot be processed
263        """
264    
265        if isinstance(arg,numarray.NumArray):
266            return arg.rank
267        elif isinstance(arg,escript.Data):
268            return arg.getRank()
269        elif isinstance(arg,float):
270            return 0
271        elif isinstance(arg,int):
272            return 0
273        elif isinstance(arg,Symbol):
274            return arg.getRank()
275        else:
276          raise TypeError,"getShape: cannot identify shape"
277    def getShape(arg):
278      """      """
279      identifies the shape of its argument      identifies the shape of its argument
280    
# Line 290  def pokeShape(arg): Line 296  def pokeShape(arg):
296      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
297          return arg.getShape()          return arg.getShape()
298      else:      else:
299        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
300    
301  def pokeDim(arg):  def pokeDim(arg):
302      """      """
# Line 313  def commonShape(arg0,arg1): Line 319  def commonShape(arg0,arg1):
319      """      """
320      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.
321    
322      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
323      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
324      @return: the shape of arg0 or arg1 such that the left port equals the shape of arg0 and the right end equals the shape of arg1.      @return: the shape of arg0 or arg1 such that the left port equals the shape of arg0 and the right end equals the shape of arg1.
325      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
326      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
327      """      """
328      sh0=pokeShape(arg0)      sh0=getShape(arg0)
329      sh1=pokeShape(arg1)      sh1=getShape(arg1)
330      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
331         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
332               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
# Line 338  def commonDim(*args): Line 344  def commonDim(*args):
344      """      """
345      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.
346    
347      @param *args: given objects      @param args: given objects
348      @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
349               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
350      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 360  def testForZero(arg): Line 366  def testForZero(arg):
366    
367      @param arg: a given object      @param arg: a given object
368      @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}
369      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
370      @rtype : C{bool}      @rtype: C{bool}
371      """      """
372      try:      if isinstance(arg,numarray.NumArray):
373           return not Lsup(arg)>0.
374        elif isinstance(arg,escript.Data):
375           return False
376        elif isinstance(arg,float):
377         return not Lsup(arg)>0.         return not Lsup(arg)>0.
378      except TypeError:      elif isinstance(arg,int):
379           return not Lsup(arg)>0.
380        elif isinstance(arg,Symbol):
381           return False
382        else:
383         return False         return False
384    
385  def matchType(arg0=0.,arg1=0.):  def matchType(arg0=0.,arg1=0.):
# Line 386  def matchType(arg0=0.,arg1=0.): Line 400  def matchType(arg0=0.,arg1=0.):
400         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
401            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
402         elif isinstance(arg1,float):         elif isinstance(arg1,float):
403            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
404         elif isinstance(arg1,int):         elif isinstance(arg1,int):
405            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
406         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
407            pass            pass
408         else:         else:
# Line 412  def matchType(arg0=0.,arg1=0.): Line 426  def matchType(arg0=0.,arg1=0.):
426         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
427            pass            pass
428         elif isinstance(arg1,float):         elif isinstance(arg1,float):
429            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
430         elif isinstance(arg1,int):         elif isinstance(arg1,int):
431            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
432         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
433            pass            pass
434         else:         else:
435            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
436      elif isinstance(arg0,float):      elif isinstance(arg0,float):
437         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
438            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
439         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
440            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
441         elif isinstance(arg1,float):         elif isinstance(arg1,float):
442            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
443            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
444         elif isinstance(arg1,int):         elif isinstance(arg1,int):
445            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
446            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
447         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
448            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
449         else:         else:
450            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
451      elif isinstance(arg0,int):      elif isinstance(arg0,int):
452         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
453            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
454         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
455            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
456         elif isinstance(arg1,float):         elif isinstance(arg1,float):
457            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
458            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
459         elif isinstance(arg1,int):         elif isinstance(arg1,int):
460            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
461            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
462         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
463            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
464         else:         else:
465            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
466      else:      else:
# Line 456  def matchType(arg0=0.,arg1=0.): Line 470  def matchType(arg0=0.,arg1=0.):
470    
471  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
472      """      """
473            return representations of arg0 amd arg1 which ahve the same shape
474    
475      If shape is not given the shape "largest" shape of args is used.      @param arg0: a given object
476        @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
477      @param args: a given ob      @param arg1: a given object
478      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
479      @return: True if the argument is identical to zero.      @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
480      @rtype: C{list} of C{int}      @rtype: C{tuple}
481      """      """
482      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
483      sh0=pokeShape(arg0)      sh0=getShape(arg0)
484      sh1=pokeShape(arg1)      sh1=getShape(arg1)
485      if len(sh0)<len(sh):      if len(sh0)<len(sh):
486         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
487      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
488         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
489      else:      else:
490         return arg0,arg1         return arg0,arg1
491  #=========================================================  #=========================================================
# Line 491  class Symbol(object): Line 505  class Symbol(object):
505         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
506         symbols or any other object.         symbols or any other object.
507    
508         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
509         @type arg: C{list}         @type args: C{list}
510         @param shape: the shape         @param shape: the shape
511         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
512         @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 549  class Symbol(object):
549         """         """
550         the shape of the symbol.         the shape of the symbol.
551    
552         @return : the shape of the symbol.         @return: the shape of the symbol.
553         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
554         """         """
555         return self.__shape         return self.__shape
# Line 544  class Symbol(object): Line 558  class Symbol(object):
558         """         """
559         the spatial dimension         the spatial dimension
560    
561         @return : the spatial dimension         @return: the spatial dimension
562         @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.
563         """         """
564         return self.__dim         return self.__dim
# Line 568  class Symbol(object): Line 582  class Symbol(object):
582         """         """
583         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.
584    
585         @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].
586         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
587         @rtype: C{list} of objects         @rtype: C{list} of objects
588         @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.
589         """         """
590         out=[]         out=[]
591         for a in self.getArgument():         for a in self.getArgument():
# Line 595  class Symbol(object): Line 609  class Symbol(object):
609            if isinstance(a,Symbol):            if isinstance(a,Symbol):
610               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
611            else:            else:
612                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
613                if len(s)>0:                if len(s)>0:
614                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
615                else:                else:
616                   out.append(a)                   out.append(a)
617         return out         return out
# Line 687  class Symbol(object): Line 701  class Symbol(object):
701         else:         else:
702            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
703            if len(s)>0:            if len(s)>0:
704               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
705            else:            else:
706               return 0.               return 0.
707    
# Line 695  class Symbol(object): Line 709  class Symbol(object):
709         """         """
710         returns -self.         returns -self.
711    
712         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
713         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
714         """         """
715         return self*(-1.)         return self*(-1.)
# Line 704  class Symbol(object): Line 718  class Symbol(object):
718         """         """
719         returns +self.         returns +self.
720    
721         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
722         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
723         """         """
724         return self*(1.)         return self*(1.)
725    
726     def __abs__(self):     def __abs__(self):
727         """         """
728         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
729         """         """
730         return Abs_Symbol(self)         return Abs_Symbol(self)
731    
# Line 721  class Symbol(object): Line 735  class Symbol(object):
735    
736         @param other: object to be added to this object         @param other: object to be added to this object
737         @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}.
738         @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}
739         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
740         """         """
741         return add(self,other)         return add(self,other)
# Line 732  class Symbol(object): Line 746  class Symbol(object):
746    
747         @param other: object this object is added to         @param other: object this object is added to
748         @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}.
749         @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
750         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
751         """         """
752         return add(other,self)         return add(other,self)
# Line 743  class Symbol(object): Line 757  class Symbol(object):
757    
758         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
759         @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}.
760         @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
761         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
762         """         """
763         return add(self,-other)         return add(self,-other)
# Line 754  class Symbol(object): Line 768  class Symbol(object):
768    
769         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
770         @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}.
771         @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}.
772         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
773         """         """
774         return add(-self,other)         return add(-self,other)
# Line 765  class Symbol(object): Line 779  class Symbol(object):
779    
780         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
781         @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}.
782         @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}.
783         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
784         """         """
785         return mult(self,other)         return mult(self,other)
# Line 776  class Symbol(object): Line 790  class Symbol(object):
790    
791         @param other: object this object is multiplied with         @param other: object this object is multiplied with
792         @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}.
793         @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.
794         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
795         """         """
796         return mult(other,self)         return mult(other,self)
# Line 787  class Symbol(object): Line 801  class Symbol(object):
801    
802         @param other: object dividing this object         @param other: object dividing this object
803         @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}.
804         @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}
805         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
806         """         """
807         return quotient(self,other)         return quotient(self,other)
# Line 798  class Symbol(object): Line 812  class Symbol(object):
812    
813         @param other: object dividing this object         @param other: object dividing this object
814         @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}.
815         @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
816         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
817         """         """
818         return quotient(other,self)         return quotient(other,self)
# Line 809  class Symbol(object): Line 823  class Symbol(object):
823    
824         @param other: exponent         @param other: exponent
825         @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}.
826         @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}
827         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
828         """         """
829         return power(self,other)         return power(self,other)
# Line 820  class Symbol(object): Line 834  class Symbol(object):
834    
835         @param other: basis         @param other: basis
836         @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}.
837         @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
838         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
839         """         """
840         return power(other,self)         return power(other,self)
841    
842       def __getitem__(self,index):
843           """
844           returns the slice defined by index
845    
846           @param index: defines a
847           @type index: C{slice} or C{int} or a C{tuple} of them
848           @return: a L{Symbol} representing the slice defined by index
849           @rtype: L{DependendSymbol}
850           """
851           return GetSlice_Symbol(self,index)
852    
853  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
854     """     """
855     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.
856     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  
857        
858     Example:     Example::
859        
860     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
861     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
862     print u1==u2       print u1==u2
863     False       False
864        
865        but       but::
866    
867     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
868     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
869     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
870     print u1==u2, u1==u3       print u1==u2, u1==u3
871     True False       True False
872    
873     @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.
874     """     """
# Line 875  class DependendSymbol(Symbol): Line 900  class DependendSymbol(Symbol):
900  #=========================================================  #=========================================================
901  #  Unary operations prserving the shape  #  Unary operations prserving the shape
902  #========================================================  #========================================================
903    class GetSlice_Symbol(DependendSymbol):
904       """
905       L{Symbol} representing getting a slice for a L{Symbol}
906       """
907       def __init__(self,arg,index):
908          """
909          initialization of wherePositive L{Symbol} with argument arg
910          @param arg: argument
911          @type arg: L{Symbol}.
912          @param index: defines index
913          @type index: C{slice} or C{int} or a C{tuple} of them
914          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
915          @raises ValueError: if a step is given
916          """
917          if not isinstance(index,tuple): index=(index,)
918          if len(index)>arg.getRank():
919               raise IndexError,"GetSlice_Symbol: index out of range."
920          sh=()
921          index2=()
922          for i in range(len(index)):
923             ix=index[i]
924             if isinstance(ix,int):
925                if ix<0 or ix>=arg.getShape()[i]:
926                   raise ValueError,"GetSlice_Symbol: index out of range."
927                index2=index2+(ix,)
928             else:
929               if not ix.step==None:
930                 raise ValueError,"GetSlice_Symbol: steping is not supported."
931               if ix.start==None:
932                  s=0
933               else:
934                  s=ix.start
935               if ix.stop==None:
936                  e=arg.getShape()[i]
937               else:
938                  e=ix.stop
939                  if e>arg.getShape()[i]:
940                     raise IndexError,"GetSlice_Symbol: index out of range."
941               index2=index2+(slice(s,e),)
942               if e>s:
943                   sh=sh+(e-s,)
944               elif s>e:
945                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
946          for i in range(len(index),arg.getRank()):
947              index2=index2+(slice(0,arg.getShape()[i]),)
948              sh=sh+(arg.getShape()[i],)
949          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
950    
951       def getMyCode(self,argstrs,format="escript"):
952          """
953          returns a program code that can be used to evaluate the symbol.
954    
955          @param argstrs: gives for each argument a string representing the argument for the evaluation.
956          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
957          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
958          @type format: C{str}
959          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
960          @rtype: C{str}
961          @raise NotImplementedError: if the requested format is not available
962          """
963          if format=="escript" or format=="str"  or format=="text":
964             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
965          else:
966             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
967    
968       def substitute(self,argvals):
969          """
970          assigns new values to symbols in the definition of the symbol.
971          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
972    
973          @param argvals: new values assigned to symbols
974          @type argvals: C{dict} with keywords of type L{Symbol}.
975          @return: result of the substitution process. Operations are executed as much as possible.
976          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
977          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
978          """
979          if argvals.has_key(self):
980             arg=argvals[self]
981             if self.isAppropriateValue(arg):
982                return arg
983             else:
984                raise TypeError,"%s: new value is not appropriate."%str(self)
985          else:
986             args=self.getSubstitutedArguments(argvals)
987             arg=args[0]
988             index=args[1]
989             return arg.__getitem__(index)
990    
991  def log10(arg):  def log10(arg):
992     """     """
993     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
994    
995     @param arg: argument     @param arg: argument
996     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
997     @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.
998     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
999     """     """
1000     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 903  def wherePositive(arg): Line 1016  def wherePositive(arg):
1016    
1017     @param arg: argument     @param arg: argument
1018     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1019     @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.
1020     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1021     """     """
1022     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1023        if arg.rank==0:        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1024           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1025             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))  
1026     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1027        return arg._wherePositive()        return arg._wherePositive()
1028     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 953  class WherePositive_Symbol(DependendSymb Line 1062  class WherePositive_Symbol(DependendSymb
1062        @type format: C{str}        @type format: C{str}
1063        @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.
1064        @rtype: C{str}        @rtype: C{str}
1065        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1066        """        """
1067        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1068            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 989  def whereNegative(arg): Line 1098  def whereNegative(arg):
1098    
1099     @param arg: argument     @param arg: argument
1100     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1101     @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.
1102     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1103     """     """
1104     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1105        if arg.rank==0:        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1106           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1107             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))  
1108     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1109        return arg._whereNegative()        return arg._whereNegative()
1110     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1039  class WhereNegative_Symbol(DependendSymb Line 1144  class WhereNegative_Symbol(DependendSymb
1144        @type format: C{str}        @type format: C{str}
1145        @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.
1146        @rtype: C{str}        @rtype: C{str}
1147        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1148        """        """
1149        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1150            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1075  def whereNonNegative(arg): Line 1180  def whereNonNegative(arg):
1180    
1181     @param arg: argument     @param arg: argument
1182     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1183     @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.
1184     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1185     """     """
1186     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1187        if arg.rank==0:        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1188           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1189             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))  
1190     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1191        return arg._whereNonNegative()        return arg._whereNonNegative()
1192     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1109  def whereNonPositive(arg): Line 1210  def whereNonPositive(arg):
1210    
1211     @param arg: argument     @param arg: argument
1212     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1213     @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.
1214     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1215     """     """
1216     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1217        if arg.rank==0:        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1218           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1219             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.  
1220     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1221        return arg._whereNonPositive()        return arg._whereNonPositive()
1222     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1145  def whereZero(arg,tol=0.): Line 1242  def whereZero(arg,tol=0.):
1242     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1243     @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.
1244     @type tol: C{float}     @type tol: C{float}
1245     @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.
1246     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1247     """     """
1248     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1249        if arg.rank==0:        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1250           if abs(arg)<=tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1251             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.  
1252     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1253        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1254     elif isinstance(arg,float):     elif isinstance(arg,float):
1255        if abs(arg)<=tol:        if abs(arg)<=tol:
1256          return 1.          return 1.
# Line 1198  class WhereZero_Symbol(DependendSymbol): Line 1288  class WhereZero_Symbol(DependendSymbol):
1288        @type format: C{str}        @type format: C{str}
1289        @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.
1290        @rtype: C{str}        @rtype: C{str}
1291        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1292        """        """
1293        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1294           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 1322  def whereNonZero(arg,tol=0.):
1322    
1323     @param arg: argument     @param arg: argument
1324     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1325     @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.
1326     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1327     """     """
1328     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1329        if arg.rank==0:        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1330          if abs(arg)>tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1331             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.  
1332     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1333        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1334     elif isinstance(arg,float):     elif isinstance(arg,float):
1335        if abs(arg)>tol:        if abs(arg)>tol:
1336          return 1.          return 1.
# Line 1263  def whereNonZero(arg,tol=0.): Line 1346  def whereNonZero(arg,tol=0.):
1346     else:     else:
1347        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1348    
1349    def erf(arg):
1350       """
1351       returns erf of argument arg
1352    
1353       @param arg: argument
1354       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1355       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1356       @raises TypeError: if the type of the argument is not expected.
1357       """
1358       if isinstance(arg,escript.Data):
1359          return arg._erf()
1360       else:
1361          raise TypeError,"erf: Unknown argument type."
1362    
1363  def sin(arg):  def sin(arg):
1364     """     """
1365     returns sine of argument arg     returns sine of argument arg
1366    
1367     @param arg: argument     @param arg: argument
1368     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1369     @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.
1370     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1371     """     """
1372     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1307  class Sin_Symbol(DependendSymbol): Line 1404  class Sin_Symbol(DependendSymbol):
1404        @type format: C{str}        @type format: C{str}
1405        @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.
1406        @rtype: C{str}        @rtype: C{str}
1407        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1408        """        """
1409        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1410            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1359  def cos(arg): Line 1456  def cos(arg):
1456    
1457     @param arg: argument     @param arg: argument
1458     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1459     @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.
1460     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1461     """     """
1462     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1397  class Cos_Symbol(DependendSymbol): Line 1494  class Cos_Symbol(DependendSymbol):
1494        @type format: C{str}        @type format: C{str}
1495        @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.
1496        @rtype: C{str}        @rtype: C{str}
1497        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1498        """        """
1499        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1500            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1449  def tan(arg): Line 1546  def tan(arg):
1546    
1547     @param arg: argument     @param arg: argument
1548     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1549     @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.
1550     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1551     """     """
1552     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1487  class Tan_Symbol(DependendSymbol): Line 1584  class Tan_Symbol(DependendSymbol):
1584        @type format: C{str}        @type format: C{str}
1585        @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.
1586        @rtype: C{str}        @rtype: C{str}
1587        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1588        """        """
1589        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1590            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1539  def asin(arg): Line 1636  def asin(arg):
1636    
1637     @param arg: argument     @param arg: argument
1638     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1639     @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.
1640     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1641     """     """
1642     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1577  class Asin_Symbol(DependendSymbol): Line 1674  class Asin_Symbol(DependendSymbol):
1674        @type format: C{str}        @type format: C{str}
1675        @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.
1676        @rtype: C{str}        @rtype: C{str}
1677        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1678        """        """
1679        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1680            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1629  def acos(arg): Line 1726  def acos(arg):
1726    
1727     @param arg: argument     @param arg: argument
1728     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1729     @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.
1730     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1731     """     """
1732     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1667  class Acos_Symbol(DependendSymbol): Line 1764  class Acos_Symbol(DependendSymbol):
1764        @type format: C{str}        @type format: C{str}
1765        @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.
1766        @rtype: C{str}        @rtype: C{str}
1767        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1768        """        """
1769        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1770            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1719  def atan(arg): Line 1816  def atan(arg):
1816    
1817     @param arg: argument     @param arg: argument
1818     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1819     @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.
1820     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1821     """     """
1822     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1757  class Atan_Symbol(DependendSymbol): Line 1854  class Atan_Symbol(DependendSymbol):
1854        @type format: C{str}        @type format: C{str}
1855        @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.
1856        @rtype: C{str}        @rtype: C{str}
1857        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1858        """        """
1859        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1860            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1809  def sinh(arg): Line 1906  def sinh(arg):
1906    
1907     @param arg: argument     @param arg: argument
1908     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1909     @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.
1910     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1911     """     """
1912     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1847  class Sinh_Symbol(DependendSymbol): Line 1944  class Sinh_Symbol(DependendSymbol):
1944        @type format: C{str}        @type format: C{str}
1945        @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.
1946        @rtype: C{str}        @rtype: C{str}
1947        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1948        """        """
1949        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1950            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1899  def cosh(arg): Line 1996  def cosh(arg):
1996    
1997     @param arg: argument     @param arg: argument
1998     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1999     @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.
2000     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2001     """     """
2002     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1937  class Cosh_Symbol(DependendSymbol): Line 2034  class Cosh_Symbol(DependendSymbol):
2034        @type format: C{str}        @type format: C{str}
2035        @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.
2036        @rtype: C{str}        @rtype: C{str}
2037        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2038        """        """
2039        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2040            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1989  def tanh(arg): Line 2086  def tanh(arg):
2086    
2087     @param arg: argument     @param arg: argument
2088     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2089     @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.
2090     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2091     """     """
2092     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2027  class Tanh_Symbol(DependendSymbol): Line 2124  class Tanh_Symbol(DependendSymbol):
2124        @type format: C{str}        @type format: C{str}
2125        @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.
2126        @rtype: C{str}        @rtype: C{str}
2127        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2128        """        """
2129        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2130            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2079  def asinh(arg): Line 2176  def asinh(arg):
2176    
2177     @param arg: argument     @param arg: argument
2178     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2179     @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.
2180     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2181     """     """
2182     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2117  class Asinh_Symbol(DependendSymbol): Line 2214  class Asinh_Symbol(DependendSymbol):
2214        @type format: C{str}        @type format: C{str}
2215        @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.
2216        @rtype: C{str}        @rtype: C{str}
2217        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2218        """        """
2219        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2220            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2169  def acosh(arg): Line 2266  def acosh(arg):
2266    
2267     @param arg: argument     @param arg: argument
2268     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2269     @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.
2270     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2271     """     """
2272     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2207  class Acosh_Symbol(DependendSymbol): Line 2304  class Acosh_Symbol(DependendSymbol):
2304        @type format: C{str}        @type format: C{str}
2305        @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.
2306        @rtype: C{str}        @rtype: C{str}
2307        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2308        """        """
2309        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2310            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2259  def atanh(arg): Line 2356  def atanh(arg):
2356    
2357     @param arg: argument     @param arg: argument
2358     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2359     @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.
2360     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2361     """     """
2362     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2297  class Atanh_Symbol(DependendSymbol): Line 2394  class Atanh_Symbol(DependendSymbol):
2394        @type format: C{str}        @type format: C{str}
2395        @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.
2396        @rtype: C{str}        @rtype: C{str}
2397        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2398        """        """
2399        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2400            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2349  def exp(arg): Line 2446  def exp(arg):
2446    
2447     @param arg: argument     @param arg: argument
2448     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2449     @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.
2450     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2451     """     """
2452     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2387  class Exp_Symbol(DependendSymbol): Line 2484  class Exp_Symbol(DependendSymbol):
2484        @type format: C{str}        @type format: C{str}
2485        @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.
2486        @rtype: C{str}        @rtype: C{str}
2487        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2488        """        """
2489        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2490            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2439  def sqrt(arg): Line 2536  def sqrt(arg):
2536    
2537     @param arg: argument     @param arg: argument
2538     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2539     @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.
2540     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2541     """     """
2542     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2477  class Sqrt_Symbol(DependendSymbol): Line 2574  class Sqrt_Symbol(DependendSymbol):
2574        @type format: C{str}        @type format: C{str}
2575        @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.
2576        @rtype: C{str}        @rtype: C{str}
2577        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2578        """        """
2579        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2580            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2529  def log(arg): Line 2626  def log(arg):
2626    
2627     @param arg: argument     @param arg: argument
2628     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2629     @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.
2630     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2631     """     """
2632     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2567  class Log_Symbol(DependendSymbol): Line 2664  class Log_Symbol(DependendSymbol):
2664        @type format: C{str}        @type format: C{str}
2665        @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.
2666        @rtype: C{str}        @rtype: C{str}
2667        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2668        """        """
2669        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2670            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2619  def sign(arg): Line 2716  def sign(arg):
2716    
2717     @param arg: argument     @param arg: argument
2718     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2719     @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.
2720     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2721     """     """
2722     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2667  class Abs_Symbol(DependendSymbol): Line 2764  class Abs_Symbol(DependendSymbol):
2764        @type format: C{str}        @type format: C{str}
2765        @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.
2766        @rtype: C{str}        @rtype: C{str}
2767        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2768        """        """
2769        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2770            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2719  def minval(arg): Line 2816  def minval(arg):
2816    
2817     @param arg: argument     @param arg: argument
2818     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2819     @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.
2820     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2821     """     """
2822     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2760  class Minval_Symbol(DependendSymbol): Line 2857  class Minval_Symbol(DependendSymbol):
2857        @type format: C{str}        @type format: C{str}
2858        @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.
2859        @rtype: C{str}        @rtype: C{str}
2860        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2861        """        """
2862        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2863            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2796  def maxval(arg): Line 2893  def maxval(arg):
2893    
2894     @param arg: argument     @param arg: argument
2895     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2896     @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.
2897     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2898     """     """
2899     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2837  class Maxval_Symbol(DependendSymbol): Line 2934  class Maxval_Symbol(DependendSymbol):
2934        @type format: C{str}        @type format: C{str}
2935        @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.
2936        @rtype: C{str}        @rtype: C{str}
2937        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2938        """        """
2939        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2940            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2873  def length(arg): Line 2970  def length(arg):
2970    
2971     @param arg: argument     @param arg: argument
2972     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2973     @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.
2974     """     """
2975     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2976    
2977    def trace(arg,axis_offset=0):
2978       """
2979       returns the trace of arg which the sum of arg[k,k] over k.
2980    
2981       @param arg: argument
2982       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2983       @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
2984                      C{axis_offset} and axis_offset+1 must be equal.
2985       @type axis_offset: C{int}
2986       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2987       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2988       """
2989       if isinstance(arg,numarray.NumArray):
2990          sh=arg.shape
2991          if len(sh)<2:
2992            raise ValueError,"rank of argument must be greater than 1"
2993          if axis_offset<0 or axis_offset>len(sh)-2:
2994            raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
2995          s1=1
2996          for i in range(axis_offset): s1*=sh[i]
2997          s2=1
2998          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2999          if not sh[axis_offset] == sh[axis_offset+1]:
3000            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3001          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
3002          out=numarray.zeros([s1,s2],numarray.Float64)
3003          for i1 in range(s1):
3004            for i2 in range(s2):
3005                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
3006          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3007          return out
3008       elif isinstance(arg,escript.Data):
3009          if arg.getRank()<2:
3010            raise ValueError,"rank of argument must be greater than 1"
3011          if axis_offset<0 or axis_offset>arg.getRank()-2:
3012            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3013          s=list(arg.getShape())        
3014          if not s[axis_offset] == s[axis_offset+1]:
3015            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3016          return arg._trace(axis_offset)
3017       elif isinstance(arg,float):
3018          raise TypeError,"illegal argument type float."
3019       elif isinstance(arg,int):
3020          raise TypeError,"illegal argument type int."
3021       elif isinstance(arg,Symbol):
3022          return Trace_Symbol(arg,axis_offset)
3023       else:
3024          raise TypeError,"Unknown argument type."
3025    
3026    class Trace_Symbol(DependendSymbol):
3027       """
3028       L{Symbol} representing the result of the trace function
3029       """
3030       def __init__(self,arg,axis_offset=0):
3031          """
3032          initialization of trace L{Symbol} with argument arg
3033          @param arg: argument of function
3034          @type arg: L{Symbol}.
3035          @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
3036                      C{axis_offset} and axis_offset+1 must be equal.
3037          @type axis_offset: C{int}
3038          """
3039          if arg.getRank()<2:
3040            raise ValueError,"rank of argument must be greater than 1"
3041          if axis_offset<0 or axis_offset>arg.getRank()-2:
3042            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3043          s=list(arg.getShape())        
3044          if not s[axis_offset] == s[axis_offset+1]:
3045            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3046          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3047    
3048       def getMyCode(self,argstrs,format="escript"):
3049          """
3050          returns a program code that can be used to evaluate the symbol.
3051    
3052          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3053          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3054          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3055          @type format: C{str}
3056          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3057          @rtype: C{str}
3058          @raise NotImplementedError: if the requested format is not available
3059          """
3060          if format=="escript" or format=="str"  or format=="text":
3061             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3062          else:
3063             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3064    
3065       def substitute(self,argvals):
3066          """
3067          assigns new values to symbols in the definition of the symbol.
3068          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3069    
3070          @param argvals: new values assigned to symbols
3071          @type argvals: C{dict} with keywords of type L{Symbol}.
3072          @return: result of the substitution process. Operations are executed as much as possible.
3073          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3074          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3075          """
3076          if argvals.has_key(self):
3077             arg=argvals[self]
3078             if self.isAppropriateValue(arg):
3079                return arg
3080             else:
3081                raise TypeError,"%s: new value is not appropriate."%str(self)
3082          else:
3083             arg=self.getSubstitutedArguments(argvals)
3084             return trace(arg[0],axis_offset=arg[1])
3085    
3086       def diff(self,arg):
3087          """
3088          differential of this object
3089    
3090          @param arg: the derivative is calculated with respect to arg
3091          @type arg: L{escript.Symbol}
3092          @return: derivative with respect to C{arg}
3093          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3094          """
3095          if arg==self:
3096             return identity(self.getShape())
3097          else:
3098             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3099    
3100    def transpose(arg,axis_offset=None):
3101       """
3102       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3103    
3104       @param arg: argument
3105       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3106       @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.
3107                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3108       @type axis_offset: C{int}
3109       @return: transpose of arg
3110       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3111       """
3112       if isinstance(arg,numarray.NumArray):
3113          if axis_offset==None: axis_offset=int(arg.rank/2)
3114          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3115       elif isinstance(arg,escript.Data):
3116          r=arg.getRank()
3117          if axis_offset==None: axis_offset=int(r/2)
3118          if axis_offset<0 or axis_offset>r:
3119            raise ValueError,"axis_offset must be between 0 and %s"%r
3120          return arg._transpose(axis_offset)
3121       elif isinstance(arg,float):
3122          if not ( axis_offset==0 or axis_offset==None):
3123            raise ValueError,"axis_offset must be 0 for float argument"
3124          return arg
3125       elif isinstance(arg,int):
3126          if not ( axis_offset==0 or axis_offset==None):
3127            raise ValueError,"axis_offset must be 0 for int argument"
3128          return float(arg)
3129       elif isinstance(arg,Symbol):
3130          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3131          return Transpose_Symbol(arg,axis_offset)
3132       else:
3133          raise TypeError,"Unknown argument type."
3134    
3135    class Transpose_Symbol(DependendSymbol):
3136       """
3137       L{Symbol} representing the result of the transpose function
3138       """
3139       def __init__(self,arg,axis_offset=None):
3140          """
3141          initialization of transpose L{Symbol} with argument arg
3142    
3143          @param arg: argument of function
3144          @type arg: L{Symbol}.
3145          @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.
3146                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3147          @type axis_offset: C{int}
3148          """
3149          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3150          if axis_offset<0 or axis_offset>arg.getRank():
3151            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3152          s=arg.getShape()
3153          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3154    
3155       def getMyCode(self,argstrs,format="escript"):
3156          """
3157          returns a program code that can be used to evaluate the symbol.
3158    
3159          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3160          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3161          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3162          @type format: C{str}
3163          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3164          @rtype: C{str}
3165          @raise NotImplementedError: if the requested format is not available
3166          """
3167          if format=="escript" or format=="str"  or format=="text":
3168             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3169          else:
3170             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3171    
3172       def substitute(self,argvals):
3173          """
3174          assigns new values to symbols in the definition of the symbol.
3175          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3176    
3177          @param argvals: new values assigned to symbols
3178          @type argvals: C{dict} with keywords of type L{Symbol}.
3179          @return: result of the substitution process. Operations are executed as much as possible.
3180          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3181          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3182          """
3183          if argvals.has_key(self):
3184             arg=argvals[self]
3185             if self.isAppropriateValue(arg):
3186                return arg
3187             else:
3188                raise TypeError,"%s: new value is not appropriate."%str(self)
3189          else:
3190             arg=self.getSubstitutedArguments(argvals)
3191             return transpose(arg[0],axis_offset=arg[1])
3192    
3193       def diff(self,arg):
3194          """
3195          differential of this object
3196    
3197          @param arg: the derivative is calculated with respect to arg
3198          @type arg: L{escript.Symbol}
3199          @return: derivative with respect to C{arg}
3200          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3201          """
3202          if arg==self:
3203             return identity(self.getShape())
3204          else:
3205             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3206    
3207    def swap_axes(arg,axis0=0,axis1=1):
3208       """
3209       returns the swap of arg by swaping the components axis0 and axis1
3210    
3211       @param arg: argument
3212       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3213       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3214       @type axis0: C{int}
3215       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3216       @type axis1: C{int}
3217       @return: C{arg} with swaped components
3218       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3219       """
3220       if axis0 > axis1:
3221          axis0,axis1=axis1,axis0
3222       if isinstance(arg,numarray.NumArray):
3223          return numarray.swapaxes(arg,axis0,axis1)
3224       elif isinstance(arg,escript.Data):
3225          return arg._swap_axes(axis0,axis1)
3226       elif isinstance(arg,float):
3227          raise TyepError,"float argument is not supported."
3228       elif isinstance(arg,int):
3229          raise TyepError,"int argument is not supported."
3230       elif isinstance(arg,Symbol):
3231          return SwapAxes_Symbol(arg,axis0,axis1)
3232       else:
3233          raise TypeError,"Unknown argument type."
3234    
3235    class SwapAxes_Symbol(DependendSymbol):
3236       """
3237       L{Symbol} representing the result of the swap function
3238       """
3239       def __init__(self,arg,axis0=0,axis1=1):
3240          """
3241          initialization of swap L{Symbol} with argument arg
3242    
3243          @param arg: argument
3244          @type arg: L{Symbol}.
3245          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3246          @type axis0: C{int}
3247          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3248          @type axis1: C{int}
3249          """
3250          if arg.getRank()<2:
3251             raise ValueError,"argument must have at least rank 2."
3252          if axis0<0 or axis0>arg.getRank()-1:
3253             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3254          if axis1<0 or axis1>arg.getRank()-1:
3255             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3256          if axis0 == axis1:
3257             raise ValueError,"axis indices must be different."
3258          if axis0 > axis1:
3259             axis0,axis1=axis1,axis0
3260          s=arg.getShape()
3261          s_out=[]
3262          for i in range(len(s)):
3263             if i == axis0:
3264                s_out.append(s[axis1])
3265             elif i == axis1:
3266                s_out.append(s[axis0])
3267             else:
3268                s_out.append(s[i])
3269          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3270    
3271       def getMyCode(self,argstrs,format="escript"):
3272          """
3273          returns a program code that can be used to evaluate the symbol.
3274    
3275          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3276          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3277          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3278          @type format: C{str}
3279          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3280          @rtype: C{str}
3281          @raise NotImplementedError: if the requested format is not available
3282          """
3283          if format=="escript" or format=="str"  or format=="text":
3284             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3285          else:
3286             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3287    
3288       def substitute(self,argvals):
3289          """
3290          assigns new values to symbols in the definition of the symbol.
3291          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3292    
3293          @param argvals: new values assigned to symbols
3294          @type argvals: C{dict} with keywords of type L{Symbol}.
3295          @return: result of the substitution process. Operations are executed as much as possible.
3296          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3297          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3298          """
3299          if argvals.has_key(self):
3300             arg=argvals[self]
3301             if self.isAppropriateValue(arg):
3302                return arg
3303             else:
3304                raise TypeError,"%s: new value is not appropriate."%str(self)
3305          else:
3306             arg=self.getSubstitutedArguments(argvals)
3307             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3308    
3309       def diff(self,arg):
3310          """
3311          differential of this object
3312    
3313          @param arg: the derivative is calculated with respect to arg
3314          @type arg: L{escript.Symbol}
3315          @return: derivative with respect to C{arg}
3316          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3317          """
3318          if arg==self:
3319             return identity(self.getShape())
3320          else:
3321             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3322    
3323    def symmetric(arg):
3324        """
3325        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3326    
3327        @param arg: square matrix. Must have rank 2 or 4 and be square.
3328        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3329        @return: symmetric part of arg
3330        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3331        """
3332        if isinstance(arg,numarray.NumArray):
3333          if arg.rank==2:
3334            if not (arg.shape[0]==arg.shape[1]):
3335               raise ValueError,"argument must be square."
3336          elif arg.rank==4:
3337            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3338               raise ValueError,"argument must be square."
3339          else:
3340            raise ValueError,"rank 2 or 4 is required."
3341          return (arg+transpose(arg))/2
3342        elif isinstance(arg,escript.Data):
3343          if arg.getRank()==2:
3344            if not (arg.getShape()[0]==arg.getShape()[1]):
3345               raise ValueError,"argument must be square."
3346            return arg._symmetric()
3347          elif arg.getRank()==4:
3348            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3349               raise ValueError,"argument must be square."
3350            return arg._symmetric()
3351          else:
3352            raise ValueError,"rank 2 or 4 is required."
3353        elif isinstance(arg,float):
3354          return arg
3355        elif isinstance(arg,int):
3356          return float(arg)
3357        elif isinstance(arg,Symbol):
3358          if arg.getRank()==2:
3359            if not (arg.getShape()[0]==arg.getShape()[1]):
3360               raise ValueError,"argument must be square."
3361          elif arg.getRank()==4:
3362            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3363               raise ValueError,"argument must be square."
3364          else:
3365            raise ValueError,"rank 2 or 4 is required."
3366          return (arg+transpose(arg))/2
3367        else:
3368          raise TypeError,"symmetric: Unknown argument type."
3369    
3370    def nonsymmetric(arg):
3371        """
3372        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3373    
3374        @param arg: square matrix. Must have rank 2 or 4 and be square.
3375        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3376        @return: nonsymmetric part of arg
3377        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3378        """
3379        if isinstance(arg,numarray.NumArray):
3380          if arg.rank==2:
3381            if not (arg.shape[0]==arg.shape[1]):
3382               raise ValueError,"nonsymmetric: argument must be square."
3383          elif arg.rank==4:
3384            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3385               raise ValueError,"nonsymmetric: argument must be square."
3386          else:
3387            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3388          return (arg-transpose(arg))/2
3389        elif isinstance(arg,escript.Data):
3390          if arg.getRank()==2:
3391            if not (arg.getShape()[0]==arg.getShape()[1]):
3392               raise ValueError,"argument must be square."
3393            return arg._nonsymmetric()
3394          elif arg.getRank()==4:
3395            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3396               raise ValueError,"argument must be square."
3397            return arg._nonsymmetric()
3398          else:
3399            raise ValueError,"rank 2 or 4 is required."
3400        elif isinstance(arg,float):
3401          return arg
3402        elif isinstance(arg,int):
3403          return float(arg)
3404        elif isinstance(arg,Symbol):
3405          if arg.getRank()==2:
3406            if not (arg.getShape()[0]==arg.getShape()[1]):
3407               raise ValueError,"nonsymmetric: argument must be square."
3408          elif arg.getRank()==4:
3409            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3410               raise ValueError,"nonsymmetric: argument must be square."
3411          else:
3412            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3413          return (arg-transpose(arg))/2
3414        else:
3415          raise TypeError,"nonsymmetric: Unknown argument type."
3416    
3417    def inverse(arg):
3418        """
3419        returns the inverse of the square matrix arg.
3420    
3421        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3422        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3423        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3424        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3425        @note: for L{escript.Data} objects the dimension is restricted to 3.
3426        """
3427        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3428        if isinstance(arg,numarray.NumArray):
3429          return numarray.linear_algebra.inverse(arg)
3430        elif isinstance(arg,escript.Data):
3431          return escript_inverse(arg)
3432        elif isinstance(arg,float):
3433          return 1./arg
3434        elif isinstance(arg,int):
3435          return 1./float(arg)
3436        elif isinstance(arg,Symbol):
3437          return Inverse_Symbol(arg)
3438        else:
3439          raise TypeError,"inverse: Unknown argument type."
3440    
3441    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3442          "arg is a Data objects!!!"
3443          if not arg.getRank()==2:
3444            raise ValueError,"escript_inverse: argument must have rank 2"
3445          s=arg.getShape()      
3446          if not s[0] == s[1]:
3447            raise ValueError,"escript_inverse: argument must be a square matrix."
3448          out=escript.Data(0.,s,arg.getFunctionSpace())
3449          if s[0]==1:
3450              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3451                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3452              out[0,0]=1./arg[0,0]
3453          elif s[0]==2:
3454              A11=arg[0,0]
3455              A12=arg[0,1]
3456              A21=arg[1,0]
3457              A22=arg[1,1]
3458              D = A11*A22-A12*A21
3459              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3460                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3461              D=1./D
3462              out[0,0]= A22*D
3463              out[1,0]=-A21*D
3464              out[0,1]=-A12*D
3465              out[1,1]= A11*D
3466          elif s[0]==3:
3467              A11=arg[0,0]
3468              A21=arg[1,0]
3469              A31=arg[2,0]
3470              A12=arg[0,1]
3471              A22=arg[1,1]
3472              A32=arg[2,1]
3473              A13=arg[0,2]
3474              A23=arg[1,2]
3475              A33=arg[2,2]
3476              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3477              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3478                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3479              D=1./D
3480              out[0,0]=(A22*A33-A23*A32)*D
3481              out[1,0]=(A31*A23-A21*A33)*D
3482              out[2,0]=(A21*A32-A31*A22)*D
3483              out[0,1]=(A13*A32-A12*A33)*D
3484              out[1,1]=(A11*A33-A31*A13)*D
3485              out[2,1]=(A12*A31-A11*A32)*D
3486              out[0,2]=(A12*A23-A13*A22)*D
3487              out[1,2]=(A13*A21-A11*A23)*D
3488              out[2,2]=(A11*A22-A12*A21)*D
3489          else:
3490             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3491          return out
3492    
3493    class Inverse_Symbol(DependendSymbol):
3494       """
3495       L{Symbol} representing the result of the inverse function
3496       """
3497       def __init__(self,arg):
3498          """
3499          initialization of inverse L{Symbol} with argument arg
3500          @param arg: argument of function
3501          @type arg: L{Symbol}.
3502          """
3503          if not arg.getRank()==2:
3504            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3505          s=arg.getShape()
3506          if not s[0] == s[1]:
3507            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3508          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3509    
3510       def getMyCode(self,argstrs,format="escript"):
3511          """
3512          returns a program code that can be used to evaluate the symbol.
3513    
3514          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3515          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3516          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3517          @type format: C{str}
3518          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3519          @rtype: C{str}
3520          @raise NotImplementedError: if the requested format is not available
3521          """
3522          if format=="escript" or format=="str"  or format=="text":
3523             return "inverse(%s)"%argstrs[0]
3524          else:
3525             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3526    
3527       def substitute(self,argvals):
3528          """
3529          assigns new values to symbols in the definition of the symbol.
3530          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3531    
3532          @param argvals: new values assigned to symbols
3533          @type argvals: C{dict} with keywords of type L{Symbol}.
3534          @return: result of the substitution process. Operations are executed as much as possible.
3535          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3536          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3537          """
3538          if argvals.has_key(self):
3539             arg=argvals[self]
3540             if self.isAppropriateValue(arg):
3541                return arg
3542             else:
3543                raise TypeError,"%s: new value is not appropriate."%str(self)
3544          else:
3545             arg=self.getSubstitutedArguments(argvals)
3546             return inverse(arg[0])
3547    
3548       def diff(self,arg):
3549          """
3550          differential of this object
3551    
3552          @param arg: the derivative is calculated with respect to arg
3553          @type arg: L{escript.Symbol}
3554          @return: derivative with respect to C{arg}
3555          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3556          """
3557          if arg==self:
3558             return identity(self.getShape())
3559          else:
3560             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3561    
3562    def eigenvalues(arg):
3563        """
3564        returns the eigenvalues of the square matrix arg.
3565    
3566        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3567                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3568        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3569        @return: the eigenvalues in increasing order.
3570        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3571        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3572        """
3573        if isinstance(arg,numarray.NumArray):
3574          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3575          out.sort()
3576          return out
3577        elif isinstance(arg,escript.Data):
3578          return arg._eigenvalues()
3579        elif isinstance(arg,Symbol):
3580          if not arg.getRank()==2:
3581            raise ValueError,"eigenvalues: argument must have rank 2"
3582          s=arg.getShape()      
3583          if not s[0] == s[1]:
3584            raise ValueError,"eigenvalues: argument must be a square matrix."
3585          if s[0]==1:
3586              return arg[0]
3587          elif s[0]==2:
3588              arg1=symmetric(arg)
3589              A11=arg1[0,0]
3590              A12=arg1[0,1]
3591              A22=arg1[1,1]
3592              trA=(A11+A22)/2.
3593              A11-=trA
3594              A22-=trA
3595              s=sqrt(A12**2-A11*A22)
3596              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3597          elif s[0]==3:
3598              arg1=symmetric(arg)
3599              A11=arg1[0,0]
3600              A12=arg1[0,1]
3601              A22=arg1[1,1]
3602              A13=arg1[0,2]
3603              A23=arg1[1,2]
3604              A33=arg1[2,2]
3605              trA=(A11+A22+A33)/3.
3606              A11-=trA
3607              A22-=trA
3608              A33-=trA
3609              A13_2=A13**2
3610              A23_2=A23**2
3611              A12_2=A12**2
3612              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3613              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3614              sq_p=sqrt(p/3.)
3615              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
3616              sq_p*=2.
3617              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3618               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3619               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3620              return trA+sq_p*f
3621          else:
3622             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3623        elif isinstance(arg,float):
3624          return arg
3625        elif isinstance(arg,int):
3626          return float(arg)
3627        else:
3628          raise TypeError,"eigenvalues: Unknown argument type."
3629    
3630    def eigenvalues_and_eigenvectors(arg):
3631        """
3632        returns the eigenvalues and eigenvectors of the square matrix arg.
3633    
3634        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3635                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3636        @type arg: L{escript.Data}
3637        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3638                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3639                 the eigenvector coresponding to the i-th eigenvalue.
3640        @rtype: L{tuple} of L{escript.Data}.
3641        @note: The dimension is restricted to 3.
3642        """
3643        if isinstance(arg,numarray.NumArray):
3644          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3645        elif isinstance(arg,escript.Data):
3646          return arg._eigenvalues_and_eigenvectors()
3647        elif isinstance(arg,Symbol):
3648          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3649        elif isinstance(arg,float):
3650          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3651        elif isinstance(arg,int):
3652          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3653        else:
3654          raise TypeError,"eigenvalues: Unknown argument type."
3655  #=======================================================  #=======================================================
3656  #  Binary operations:  #  Binary operations:
3657  #=======================================================  #=======================================================
# Line 2920  class Add_Symbol(DependendSymbol): Line 3695  class Add_Symbol(DependendSymbol):
3695         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3696         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3697         """         """
3698         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3699         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3700         if not sh0==sh1:         if not sh0==sh1:
3701            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3702         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 2936  class Add_Symbol(DependendSymbol): Line 3711  class Add_Symbol(DependendSymbol):
3711        @type format: C{str}        @type format: C{str}
3712        @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.
3713        @rtype: C{str}        @rtype: C{str}
3714        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3715        """        """
3716        if format=="str" or format=="text":        if format=="str" or format=="text":
3717           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 2995  def mult(arg0,arg1): Line 3770  def mult(arg0,arg1):
3770         """         """
3771         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3772         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3773            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3774         else:         else:
3775            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3776                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3019  class Mult_Symbol(DependendSymbol): Line 3794  class Mult_Symbol(DependendSymbol):
3794         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3795         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3796         """         """
3797         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3798         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3799         if not sh0==sh1:         if not sh0==sh1:
3800            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3801         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3035  class Mult_Symbol(DependendSymbol): Line 3810  class Mult_Symbol(DependendSymbol):
3810        @type format: C{str}        @type format: C{str}
3811        @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.
3812        @rtype: C{str}        @rtype: C{str}
3813        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3814        """        """
3815        if format=="str" or format=="text":        if format=="str" or format=="text":
3816           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3095  def quotient(arg0,arg1): Line 3870  def quotient(arg0,arg1):
3870         """         """
3871         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3872         if testForZero(args[0]):         if testForZero(args[0]):
3873            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3874         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3875            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3876               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3124  class Quotient_Symbol(DependendSymbol): Line 3899  class Quotient_Symbol(DependendSymbol):
3899         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3900         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3901         """         """
3902         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3903         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3904         if not sh0==sh1:         if not sh0==sh1:
3905            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3906         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3140  class Quotient_Symbol(DependendSymbol): Line 3915  class Quotient_Symbol(DependendSymbol):
3915        @type format: C{str}        @type format: C{str}
3916        @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.
3917        @rtype: C{str}        @rtype: C{str}
3918        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3919        """        """
3920        if format=="str" or format=="text":        if format=="str" or format=="text":
3921           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3201  def power(arg0,arg1): Line 3976  def power(arg0,arg1):
3976         """         """
3977         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3978         if testForZero(args[0]):         if testForZero(args[0]):
3979            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3980         elif testForZero(args[1]):         elif testForZero(args[1]):
3981            return numarray.ones(args[0],numarray.Float)            return numarray.ones(getShape(args[1]),numarray.Float64)
3982         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3983            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3984         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 3226  class Power_Symbol(DependendSymbol): Line 4001  class Power_Symbol(DependendSymbol):
4001         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4002         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4003         """         """
4004         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4005         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4006         if not sh0==sh1:         if not sh0==sh1:
4007            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4008         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3244  class Power_Symbol(DependendSymbol): Line 4019  class Power_Symbol(DependendSymbol):
4019        @type format: C{str}        @type format: C{str}
4020        @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.
4021        @rtype: C{str}        @rtype: C{str}
4022        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4023        """        """
4024        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4025           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3304  def maximum(*args): Line 4079  def maximum(*args):
4079         if out==None:         if out==None:
4080            out=a            out=a
4081         else:         else:
4082            m=whereNegative(out-a)            diff=add(a,-out)
4083            out=m*a+(1.-m)*out            out=add(out,mult(wherePositive(diff),diff))
4084      return out      return out
4085        
4086  def minimum(*arg):  def minimum(*args):
4087      """      """
4088      the minimum over arguments args      the minimum over arguments args
4089    
# Line 3322  def minimum(*arg): Line 4097  def minimum(*arg):
4097         if out==None:         if out==None:
4098            out=a            out=a
4099         else:         else:
4100            m=whereNegative(out-a)            diff=add(a,-out)
4101            out=m*out+(1.-m)*a            out=add(out,mult(whereNegative(diff),diff))
4102      return out      return out
4103    
4104    def clip(arg,minval=None,maxval=None):
4105        """
4106        cuts the values of arg between minval and maxval
4107    
4108        @param arg: argument
4109        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4110        @param minval: lower range. If None no lower range is applied
4111        @type minval: C{float} or C{None}
4112        @param maxval: upper range. If None no upper range is applied
4113        @type maxval: C{float} or C{None}
4114        @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4115                 less then maxval are unchanged.
4116        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4117        @raise ValueError: if minval>maxval
4118        """
4119        if not minval==None and not maxval==None:
4120           if minval>maxval:
4121              raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4122        if minval == None:
4123            tmp=arg
4124        else:
4125            tmp=maximum(minval,arg)
4126        if maxval == None:
4127            return tmp
4128        else:
4129            return minimum(tmp,maxval)
4130    
4131        
4132  def inner(arg0,arg1):  def inner(arg0,arg1):
4133      """      """
# Line 3340  def inner(arg0,arg1): Line 4143  def inner(arg0,arg1):
4143      @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}
4144      @param arg1: second argument      @param arg1: second argument
4145      @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}
4146      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4147      @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
4148      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4149      """      """
4150      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4151      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4152      if not sh0==sh1:      if not sh0==sh1:
4153          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4154      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4155    
4156    def outer(arg0,arg1):
4157        """
4158        the outer product of the two argument:
4159    
4160        out[t,s]=arg0[t]*arg1[s]
4161    
4162        where
4163    
4164            - s runs through arg0.Shape
4165            - t runs through arg1.Shape
4166    
4167        @param arg0: first argument
4168        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4169        @param arg1: second argument
4170        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4171        @return: the outer product of arg0 and arg1 at each data point
4172        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4173        """
4174        return generalTensorProduct(arg0,arg1,axis_offset=0)
4175    
4176  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4177      """      """
4178        see L{matrix_mult}
4179        """
4180        return matrix_mult(arg0,arg1)
4181    
4182    def matrix_mult(arg0,arg1):
4183        """
4184      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4185    
4186      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4187    
4188            or      or
4189    
4190      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4191    
4192      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.
4193    
4194      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4195      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3370  def matrixmult(arg0,arg1): Line 4199  def matrixmult(arg0,arg1):
4199      @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
4200      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4201      """      """
4202      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4203      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4204      if not len(sh0)==2 :      if not len(sh0)==2 :
4205          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4206      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4207          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4208      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4209    
4210  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4211      """      """
4212      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  
4213      """      """
4214      return generalTensorProduct(arg0,arg1,offset=0)      return tensor_mult(arg0,arg1)
4215    
4216    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4217      """      """
4218      the tensor product of the two argument:      the tensor product of the two argument:
   
4219            
4220      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4221    
4222      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4223    
4224                   or      or
4225    
4226      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4227    
# Line 3415  def tensormult(arg0,arg1): Line 4230  def tensormult(arg0,arg1):
4230    
4231      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]
4232                                
4233                   or      or
4234    
4235      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]
4236    
4237                   or      or
4238    
4239      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]
4240    
4241      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  
4242      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.
4243    
4244      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4245      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3433  def tensormult(arg0,arg1): Line 4248  def tensormult(arg0,arg1):
4248      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4249      @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
4250      """      """
4251      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4252      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4253      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4254         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4255      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):
4256         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4257      else:      else:
4258          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4259    
4260  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4261      """      """
4262      generalized tensor product      generalized tensor product
4263    
4264      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4265    
4266      where s runs through arg0.Shape[:arg0.Rank-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[offset:]  
4267    
4268      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]
4269      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4270            - t runs through arg1.Shape[axis_offset:]
4271    
4272      @param arg0: first argument      @param arg0: first argument
4273      @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}
4274      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4275      @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}
4276      @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.
4277      @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 4281  def generalTensorProduct(arg0,arg1,offse
4281      # 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
4282      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4283         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4284             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4285         else:         else:
4286             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4287                 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)
4288             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4289             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4290             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4291             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4292             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4293             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4294             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4295             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4296             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4297             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4298             for i0 in range(d0):             for i0 in range(d0):
4299                      for i1 in range(d1):                      for i1 in range(d1):
4300                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4301             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4302             return out             return out
4303      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4304         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4305             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4306         else:         else:
4307             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)
4308      else:            else:      
4309         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4310                                    
4311  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4312     """     """
4313     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4314     """     """
4315     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4316         """         """
4317         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4318    
4319         @param arg0: numerator         @param arg0: first argument
4320         @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}.
4321         @param arg1: denominator         @param arg1: second argument
4322         @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}.
4323         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4324         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4325         """         """
4326         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4327         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4328         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4329         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4330         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4331         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4332         if not sh01==sh10:         if not sh01==sh10:
4333             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)
4334         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])
4335    
4336     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4337        """        """
# Line 3529  class GeneralTensorProduct_Symbol(Depend Line 4343  class GeneralTensorProduct_Symbol(Depend
4343        @type format: C{str}        @type format: C{str}
4344        @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.
4345        @rtype: C{str}        @rtype: C{str}
4346        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4347        """        """
4348        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4349           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])
4350        else:        else:
4351           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)
4352    
# Line 3557  class GeneralTensorProduct_Symbol(Depend Line 4371  class GeneralTensorProduct_Symbol(Depend
4371           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4372           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4373    
4374  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4375      "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!!!"
4376      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4377      shape0=arg0.getShape()[:arg0.getRank()-offset]  
4378      shape01=arg0.getShape()[arg0.getRank()-offset:]  def transposed_matrix_mult(arg0,arg1):
4379      shape10=arg1.getShape()[:offset]      """
4380      shape1=arg1.getShape()[offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4381      if not shape01==shape10:  
4382          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]
4383    
4384      # create return value:      or
4385      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())  
4386      #      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4387      s0=[[]]  
4388      for k in shape0:      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4389            s=[]  
4390            for j in s0:      The first dimension of arg0 and arg1 must match.
4391                  for i in range(k): s.append(j+[slice(i,i)])  
4392            s0=s      @param arg0: first argument of rank 2
4393      s1=[[]]      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4394      for k in shape1:      @param arg1: second argument of at least rank 1
4395            s=[]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4396            for j in s1:      @return: the product of the transposed of arg0 and arg1 at each data point
4397                  for i in range(k): s.append(j+[slice(i,i)])      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4398            s1=s      @raise ValueError: if the shapes of the arguments are not appropriate
4399      s01=[[]]      """
4400      for k in shape01:      sh0=getShape(arg0)
4401            s=[]      sh1=getShape(arg1)
4402            for j in s01:      if not len(sh0)==2 :
4403                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"first argument must have rank 2"
4404            s01=s      if not len(sh1)==2 and not len(sh1)==1:
4405            raise ValueError,"second argument must have rank 1 or 2"
4406      for i0 in s0:      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4407         for i1 in s1:  
4408           s=escript.Scalar(0.,arg0.getFunctionSpace())  def transposed_tensor_mult(arg0,arg1):
4409           for i01 in s01:      """
4410              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      the tensor product of the transposed of the first and the second argument
4411           out.__setitem__(tuple(i0+i1),s)      
4412      return out      for arg0 of rank 2 this is
4413    
4414        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4415    
4416        or
4417    
4418        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4419    
4420      
4421        and for arg0 of rank 4 this is
4422    
4423        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4424                  
4425        or
4426    
4427        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4428    
4429        or
4430    
4431        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4432    
4433        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4434        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4435    
4436        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4437    
4438        @param arg0: first argument of rank 2 or 4
4439        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4440        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4441        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4442        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4443        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4444        """
4445        sh0=getShape(arg0)
4446        sh1=getShape(arg1)
4447        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4448           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4449        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4450           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4451        else:
4452            raise ValueError,"first argument must have rank 2 or 4"
4453    
4454    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4455        """
4456        generalized tensor product of transposed of arg0 and arg1:
4457    
4458        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4459    
4460        where
4461    
4462            - s runs through arg0.Shape[axis_offset:]
4463            - r runs trough arg0.Shape[:axis_offset]
4464            - t runs through arg1.Shape[axis_offset:]
4465    
4466        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4467        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4468    
4469        @param arg0: first argument
4470        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4471        @param arg1: second argument
4472        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4473        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4474        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4475        """
4476        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4477        arg0,arg1=matchType(arg0,arg1)
4478        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4479        if isinstance(arg0,numarray.NumArray):
4480           if isinstance(arg1,Symbol):
4481               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4482           else:
4483               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4484                   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)
4485               arg0_c=arg0.copy()
4486               arg1_c=arg1.copy()
4487               sh0,sh1=arg0.shape,arg1.shape
4488               d0,d1,d01=1,1,1
4489               for i in sh0[axis_offset:]: d0*=i
4490               for i in sh1[axis_offset:]: d1*=i
4491               for i in sh0[:axis_offset]: d01*=i
4492               arg0_c.resize((d01,d0))
4493               arg1_c.resize((d01,d1))
4494               out=numarray.zeros((d0,d1),numarray.Float64)
4495               for i0 in range(d0):
4496                        for i1 in range(d1):
4497                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4498               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4499               return out
4500        elif isinstance(arg0,escript.Data):
4501           if isinstance(arg1,Symbol):
4502               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4503           else:
4504               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4505        else:      
4506           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4507                    
4508    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4509       """
4510       Symbol representing the general tensor product of the transposed of arg0 and arg1
4511       """
4512       def __init__(self,arg0,arg1,axis_offset=0):
4513           """
4514           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4515    
4516           @param arg0: first argument
4517           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4518           @param arg1: second argument
4519           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4520           @raise ValueError: inconsistent dimensions of arguments.
4521           @note: if both arguments have a spatial dimension, they must equal.
4522           """
4523           sh_arg0=getShape(arg0)
4524           sh_arg1=getShape(arg1)
4525           sh01=sh_arg0[:axis_offset]
4526           sh10=sh_arg1[:axis_offset]
4527           sh0=sh_arg0[axis_offset:]
4528           sh1=sh_arg1[axis_offset:]
4529           if not sh01==sh10:
4530               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)
4531           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4532    
4533       def getMyCode(self,argstrs,format="escript"):
4534          """
4535          returns a program code that can be used to evaluate the symbol.
4536    
4537          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4538          @type argstrs: C{list} of length 2 of C{str}.
4539          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4540          @type format: C{str}
4541          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4542          @rtype: C{str}
4543          @raise NotImplementedError: if the requested format is not available
4544          """
4545          if format=="escript" or format=="str" or format=="text":
4546             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4547          else:
4548             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4549    
4550       def substitute(self,argvals):
4551          """
4552          assigns new values to symbols in the definition of the symbol.
4553          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4554    
4555          @param argvals: new values assigned to symbols
4556          @type argvals: C{dict} with keywords of type L{Symbol}.
4557          @return: result of the substitution process. Operations are executed as much as possible.
4558          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4559          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4560          """
4561          if argvals.has_key(self):
4562             arg=argvals[self]
4563             if self.isAppropriateValue(arg):
4564                return arg
4565             else:
4566                raise TypeError,"%s: new value is not appropriate."%str(self)
4567          else:
4568             args=self.getSubstitutedArguments(argvals)
4569             return generalTransposedTensorProduct(args[0],args[1],args[2])
4570    
4571    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4572        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4573        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4574    
4575    def matrix_transposed_mult(arg0,arg1):
4576        """
4577        matrix-transposed(matrix) product of the two argument:
4578    
4579        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4580    
4581        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4582    
4583        The last dimensions of arg0 and arg1 must match.
4584    
4585        @param arg0: first argument of rank 2
4586        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4587        @param arg1: second argument of rank 2
4588        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4589        @return: the product of arg0 and the transposed of arg1 at each data point
4590        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4591        @raise ValueError: if the shapes of the arguments are not appropriate
4592        """
4593        sh0=getShape(arg0)
4594        sh1=getShape(arg1)
4595        if not len(sh0)==2 :
4596            raise ValueError,"first argument must have rank 2"
4597        if not len(sh1)==2 and not len(sh1)==1:
4598            raise ValueError,"second argument must have rank 1 or 2"
4599        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4600    
4601    def tensor_transposed_mult(arg0,arg1):
4602        """
4603        the tensor product of the first and the transpose of the second argument
4604        
4605        for arg0 of rank 2 this is
4606    
4607        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4608    
4609        and for arg0 of rank 4 this is
4610    
4611        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4612                  
4613        or
4614    
4615        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4616    
4617        In the first case the the second dimension of arg0 and arg1 must match and  
4618        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4619    
4620        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4621    
4622        @param arg0: first argument of rank 2 or 4
4623        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4624        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4625        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4626        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4627        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4628        """
4629        sh0=getShape(arg0)
4630        sh1=getShape(arg1)
4631        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4632           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4633        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4634           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4635        else:
4636            raise ValueError,"first argument must have rank 2 or 4"
4637    
4638    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4639        """
4640        generalized tensor product of transposed of arg0 and arg1:
4641    
4642        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4643    
4644        where
4645    
4646            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4647            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4648            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4649    
4650        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4651        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4652    
4653        @param arg0: first argument
4654        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4655        @param arg1: second argument
4656        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4657        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4658        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4659        """
4660        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4661        arg0,arg1=matchType(arg0,arg1)
4662        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4663        if isinstance(arg0,numarray.NumArray):
4664           if isinstance(arg1,Symbol):
4665               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4666           else:
4667               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4668                   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)
4669               arg0_c=arg0.copy()
4670               arg1_c=arg1.copy()
4671               sh0,sh1=arg0.shape,arg1.shape
4672               d0,d1,d01=1,1,1
4673               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4674               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4675               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4676               arg0_c.resize((d0,d01))
4677               arg1_c.resize((d1,d01))
4678               out=numarray.zeros((d0,d1),numarray.Float64)
4679               for i0 in range(d0):
4680                        for i1 in range(d1):
4681                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4682               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4683               return out
4684        elif isinstance(arg0,escript.Data):
4685           if isinstance(arg1,Symbol):
4686               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4687           else:
4688               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4689        else:      
4690           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4691                    
4692    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4693       """
4694       Symbol representing the general tensor product of arg0 and the transpose of arg1
4695       """
4696       def __init__(self,arg0,arg1,axis_offset=0):
4697           """
4698           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4699    
4700           @param arg0: first argument
4701           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4702           @param arg1: second argument
4703           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4704           @raise ValueError: inconsistent dimensions of arguments.
4705           @note: if both arguments have a spatial dimension, they must equal.
4706           """
4707           sh_arg0=getShape(arg0)
4708           sh_arg1=getShape(arg1)
4709           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4710           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4711           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4712           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4713           if not sh01==sh10:
4714               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)
4715           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4716    
4717       def getMyCode(self,argstrs,format="escript"):
4718          """
4719          returns a program code that can be used to evaluate the symbol.
4720    
4721          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4722          @type argstrs: C{list} of length 2 of C{str}.
4723          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4724          @type format: C{str}
4725          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4726          @rtype: C{str}
4727          @raise NotImplementedError: if the requested format is not available
4728          """
4729          if format=="escript" or format=="str" or format=="text":
4730             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4731          else:
4732             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4733    
4734       def substitute(self,argvals):
4735          """
4736          assigns new values to symbols in the definition of the symbol.
4737          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4738    
4739          @param argvals: new values assigned to symbols
4740          @type argvals: C{dict} with keywords of type L{Symbol}.
4741          @return: result of the substitution process. Operations are executed as much as possible.
4742          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4743          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4744          """
4745          if argvals.has_key(self):
4746             arg=argvals[self]
4747             if self.isAppropriateValue(arg):
4748                return arg
4749             else:
4750                raise TypeError,"%s: new value is not appropriate."%str(self)
4751          else:
4752             args=self.getSubstitutedArguments(argvals)
4753             return generalTensorTransposedProduct(args[0],args[1],args[2])
4754    
4755    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4756        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4757        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4758    
4759  #=========================================================  #=========================================================
4760  #   some little helpers  #  functions dealing with spatial dependency
4761  #=========================================================  #=========================================================
4762  def grad(arg,where=None):  def grad(arg,where=None):
4763      """      """
4764      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4765    
4766        If C{g} is the returned object, then
4767    
4768      @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.
4769                    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.
4770          - 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.
4771          - 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.
4772    
4773        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4774        @type arg: L{escript.Data} or L{Symbol}
4775      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4776                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4777        @type where: C{None} or L{escript.FunctionSpace}
4778        @return: gradient of arg.
4779        @rtype: L{escript.Data} or L{Symbol}
4780      """      """
4781      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4782         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3617  def grad(arg,where=None): Line 4786  def grad(arg,where=None):
4786         else:         else:
4787            return arg._grad(where)            return arg._grad(where)
4788      else:      else:
4789        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4790    
4791    class Grad_Symbol(DependendSymbol):
4792       """
4793       L{Symbol} representing the result of the gradient operator
4794       """
4795       def __init__(self,arg,where=None):
4796          """
4797          initialization of gradient L{Symbol} with argument arg
4798          @param arg: argument of function
4799          @type arg: L{Symbol}.
4800          @param where: FunctionSpace in which the gradient will be calculated.
4801                      If not present or C{None} an appropriate default is used.
4802          @type where: C{None} or L{escript.FunctionSpace}
4803          """
4804          d=arg.getDim()
4805          if d==None:
4806             raise ValueError,"argument must have a spatial dimension"
4807          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4808    
4809       def getMyCode(self,argstrs,format="escript"):
4810          """
4811          returns a program code that can be used to evaluate the symbol.
4812    
4813          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4814          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4815          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4816          @type format: C{str}
4817          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4818          @rtype: C{str}
4819          @raise NotImplementedError: if the requested format is not available
4820          """
4821          if format=="escript" or format=="str"  or format=="text":
4822             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4823          else:
4824             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4825    
4826       def substitute(self,argvals):
4827          """
4828          assigns new values to symbols in the definition of the symbol.
4829          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4830    
4831          @param argvals: new values assigned to symbols
4832          @type argvals: C{dict} with keywords of type L{Symbol}.
4833          @return: result of the substitution process. Operations are executed as much as possible.
4834          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4835          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4836          """
4837          if argvals.has_key(self):
4838             arg=argvals[self]
4839             if self.isAppropriateValue(arg):
4840                return arg
4841             else:
4842                raise TypeError,"%s: new value is not appropriate."%str(self)
4843          else:
4844             arg=self.getSubstitutedArguments(argvals)
4845             return grad(arg[0],where=arg[1])
4846    
4847       def diff(self,arg):
4848          """
4849          differential of this object
4850    
4851          @param arg: the derivative is calculated with respect to arg
4852          @type arg: L{escript.Symbol}
4853          @return: derivative with respect to C{arg}
4854          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4855          """
4856          if arg==self:
4857             return identity(self.getShape())
4858          else:
4859             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4860    
4861  def integrate(arg,where=None):  def integrate(arg,where=None):
4862      """      """
4863      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}
4864      its domain.      before integration.
4865    
4866      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4867        @type arg: L{escript.Data} or L{Symbol}
4868      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4869                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4870        @type where: C{None} or L{escript.FunctionSpace}
4871        @return: integral of arg.
4872        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4873      """      """
4874      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):  
4875         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4876      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4877         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 4882  def integrate(arg,where=None):
4882      else:      else:
4883        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4884    
4885    class Integrate_Symbol(DependendSymbol):
4886       """
4887       L{Symbol} representing the result of the spatial integration operator
4888       """
4889       def __init__(self,arg,where=None):
4890          """
4891          initialization of integration L{Symbol} with argument arg
4892          @param arg: argument of the integration
4893          @type arg: L{Symbol}.
4894          @param where: FunctionSpace in which the integration will be calculated.
4895                      If not present or C{None} an appropriate default is used.
4896          @type where: C{None} or L{escript.FunctionSpace}
4897          """
4898          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4899    
4900       def getMyCode(self,argstrs,format="escript"):
4901          """
4902          returns a program code that can be used to evaluate the symbol.
4903    
4904          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4905          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4906          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4907          @type format: C{str}
4908          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4909          @rtype: C{str}
4910          @raise NotImplementedError: if the requested format is not available
4911          """
4912          if format=="escript" or format=="str"  or format=="text":
4913             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4914          else:
4915             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4916    
4917       def substitute(self,argvals):
4918          """
4919          assigns new values to symbols in the definition of the symbol.
4920          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4921    
4922          @param argvals: new values assigned to symbols
4923          @type argvals: C{dict} with keywords of type L{Symbol}.
4924          @return: result of the substitution process. Operations are executed as much as possible.
4925          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4926          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4927          """
4928          if argvals.has_key(self):
4929             arg=argvals[self]
4930             if self.isAppropriateValue(arg):
4931                return arg
4932             else:
4933                raise TypeError,"%s: new value is not appropriate."%str(self)
4934          else:
4935             arg=self.getSubstitutedArguments(argvals)
4936             return integrate(arg[0],where=arg[1])
4937    
4938       def diff(self,arg):
4939          """
4940          differential of this object
4941    
4942          @param arg: the derivative is calculated with respect to arg
4943          @type arg: L{escript.Symbol}
4944          @return: derivative with respect to C{arg}
4945          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4946          """
4947          if arg==self:
4948             return identity(self.getShape())
4949          else:
4950             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4951    
4952    
4953  def interpolate(arg,where):  def interpolate(arg,where):
4954      """      """
4955      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
4956    
4957      @param arg:    interpolant      @param arg: interpolant
4958      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
4959        @param where: FunctionSpace to be interpolated to
4960        @type where: L{escript.FunctionSpace}
4961        @return: interpolated argument
4962        @rtype: C{escript.Data} or L{Symbol}
4963      """      """
4964      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4965         return Interpolated_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
4966      else:      else:
4967         return escript.Data(arg,where)         return escript.Data(arg,where)
4968    
4969  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
4970      """     """
4971      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
4972       """
4973       def __init__(self,arg,where):
4974          """
4975          initialization of interpolation L{Symbol} with argument arg
4976          @param arg: argument of the interpolation
4977          @type arg: L{Symbol}.
4978          @param where: FunctionSpace into which the argument is interpolated.
4979          @type where: L{escript.FunctionSpace}
4980          """
4981          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4982    
4983      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
4984                    be calculated.        """
4985      @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)  
4986    
4987  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
4988      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4989      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.
4990          @type format: C{str}
4991          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4992          @rtype: C{str}
4993          @raise NotImplementedError: if the requested format is not available
4994          """
4995          if format=="escript" or format=="str"  or format=="text":
4996             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
4997          else:
4998             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4999    
5000      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
5001                    to be calculated.        """
5002      """        assigns new values to symbols in the definition of the symbol.
5003      d=arg.getDomain()        The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
     return arg.interpolate(escript.FunctionOnContactOne(d))-arg.interpolate(escript.FunctionOnContactZero(d))  
5004    
5005  #=============================        @param argvals: new values assigned to symbols
5006  #        @type argvals: C{dict} with keywords of type L{Symbol}.
5007  # 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.
5008  # 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
5009  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
5010          """
5011          if argvals.has_key(self):
5012             arg=argvals[self]
5013             if self.isAppropriateValue(arg):
5014                return arg
5015             else:
5016                raise TypeError,"%s: new value is not appropriate."%str(self)
5017          else:
5018             arg=self.getSubstitutedArguments(argvals)
5019             return interpolate(arg[0],where=arg[1])
5020    
5021       def diff(self,arg):
5022          """
5023          differential of this object
5024    
5025          @param arg: the derivative is calculated with respect to arg
5026          @type arg: L{escript.Symbol}
5027          @return: derivative with respect to C{arg}
5028          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
5029          """
5030          if arg==self:
5031             return identity(self.getShape())
5032          else:
5033             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5034    
 # functions involving the underlying Domain:  
5035    
5036  def transpose(arg,axis=None):  def div(arg,where=None):
5037      """      """
5038      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
5039    
5040      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
5041        @type arg: L{escript.Data} or L{Symbol}
5042        @param where: FunctionSpace in which the divergence will be calculated.
5043                      If not present or C{None} an appropriate default is used.
5044        @type where: C{None} or L{escript.FunctionSpace}
5045        @return: divergence of arg.
5046        @rtype: L{escript.Data} or L{Symbol}
5047      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
5048      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5049         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
5050      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
5051         # 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)  
5052      else:      else:
5053         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
5054        if not arg.getShape()==(dim,):
5055          raise ValueError,"div: expected shape is (%s,)"%dim
5056        return trace(grad(arg,where))
5057    
5058  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
5059      """      """
5060      Return      returns the jump of arg across the continuity of the domain
5061    
5062      @param arg:      @param arg: argument
5063        @type arg: L{escript.Data} or L{Symbol}
5064        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5065                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5066        @type domain: C{None} or L{escript.Domain}
5067        @return: jump of arg
5068        @rtype: L{escript.Data} or L{Symbol}
5069      """      """
5070      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
5071         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)  
5072    
5073    def L2(arg):
5074        """
5075        returns the L2 norm of arg at where
5076        
5077        @param arg: function which L2 to be calculated.
5078        @type arg: L{escript.Data} or L{Symbol}
5079        @return: L2 norm of arg.
5080        @rtype: L{float} or L{Symbol}
5081        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5082        """
5083        return sqrt(integrate(inner(arg,arg)))
5084    #=============================
5085    #
5086    
5087  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5088      """      """
5089      resorts the component of arg according to index      resorts the component of arg according to index
5090    
5091      """      """
5092      pass      raise NotImplementedError
5093  #  #
5094  # $Log: util.py,v $  # $Log: util.py,v $
5095  # 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.353  
changed lines
  Added in v.1042

  ViewVC Help
Powered by ViewVC 1.1.26