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

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

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

revision 396 by gross, Wed Dec 21 05:08:25 2005 UTC revision 1044 by gross, Mon Mar 19 07:29:31 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 clip(arg,minval,maxval)  
   
 # def transpose(arg,axis=None):  
 # def trace(arg,axis0=0,axis1=1):  
 # def reorderComponents(arg,index):  
   
 # def integrate(arg,where=None):  
 # def interpolate(arg,where):  
 # def div(arg,where=None):  
 # def grad(arg,where=None):  
   
 #  
 # slicing: get  
 #          set  
 #  
 # and derivatives  
30    
31  #=========================================================  #=========================================================
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: L{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 insertTagNames(domain,**kwargs):
47        """
48        inserts tag names into the domain
49    
50        @param domain: a domain object
51        @type domain: C{escript.Domain}
52        @keyword <tag name>: tag key assigned to <tag name>
53        @type <tag name>: C{int}
54        """
55        for  k in kwargs:
56             domain.setTagMap(k,kwargs[k])
57    
58    def insertTaggedValues(target,**kwargs):
59        """
60        inserts tagged values into the tagged using tag names
61    
62        @param target: data to be filled by tagged values
63        @type target: L{escript.Data}
64        @keyword <tag name>: value to be used for <tag name>
65        @type <tag name>: C{float} or {numarray.NumArray}
66        @return: C{target}
67        @rtype: L{escript.Data}
68        """
69        for k in kwargs:
70            target.setTaggedValue(k,kwargs[k])
71        return target
72    
73        
74  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
75      """      """
76      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.
77    
78      Example:      Example::
79    
80         tmp=Scalar(..)         tmp=Scalar(..)
81         v=Vector(..)         v=Vector(..)
# Line 97  def saveDX(filename,domain=None,**data): Line 103  def saveDX(filename,domain=None,**data):
103      """      """
104      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.
105    
106      Example:      Example::
107    
108         tmp=Scalar(..)         tmp=Scalar(..)
109         v=Vector(..)         v=Vector(..)
# Line 126  def kronecker(d=3): Line 132  def kronecker(d=3):
132     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
133    
134     @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
135     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
136     @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
137     @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}  
138     """     """
139     return identityTensor(d)     return identityTensor(d)
140    
# Line 144  def identity(shape=()): Line 149  def identity(shape=()):
149     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
150     """     """
151     if len(shape)>0:     if len(shape)>0:
152        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
153        if len(shape)==1:        if len(shape)==1:
154            for i0 in range(shape[0]):            for i0 in range(shape[0]):
155               out[i0,i0]=1.               out[i0,i0]=1.
   
156        elif len(shape)==2:        elif len(shape)==2:
157            for i0 in range(shape[0]):            for i0 in range(shape[0]):
158               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 164  def identityTensor(d=3): Line 168  def identityTensor(d=3):
168     return the dxd identity matrix     return the dxd identity matrix
169    
170     @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
171     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
172     @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
173     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
174     """     """
175     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
176        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
177     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
178           return identity((d.getDim(),))
179       else:
180           return identity((d,))
181    
182  def identityTensor4(d=3):  def identityTensor4(d=3):
183     """     """
# Line 179  def identityTensor4(d=3): Line 186  def identityTensor4(d=3):
186     @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
187     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
188     @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
189     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
190     """     """
191     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
192        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
193     return identity((d,d))     elif isinstance(d,escript.Domain):
194           return identity((d.getDim(),d.getDim()))
195       else:
196           return identity((d,d))
197    
198  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
199     """     """
# Line 192  def unitVector(i=0,d=3): Line 202  def unitVector(i=0,d=3):
202     @param i: index     @param i: index
203     @type i: C{int}     @type i: C{int}
204     @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
205     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
206     @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
207     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
208     """     """
209     return kronecker(d)[i]     return kronecker(d)[i]
210    
# Line 250  def inf(arg): Line 260  def inf(arg):
260    
261      @param arg: argument      @param arg: argument
262      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
263      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
264      @rtype: C{float}      @rtype: C{float}
265      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
266      """      """
# Line 269  def inf(arg): Line 279  def inf(arg):
279  #=========================================================================  #=========================================================================
280  #   some little helpers  #   some little helpers
281  #=========================================================================  #=========================================================================
282  def pokeShape(arg):  def getRank(arg):
283        """
284        identifies the rank of its argument
285    
286        @param arg: a given object
287        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
288        @return: the rank of the argument
289        @rtype: C{int}
290        @raise TypeError: if type of arg cannot be processed
291        """
292    
293        if isinstance(arg,numarray.NumArray):
294            return arg.rank
295        elif isinstance(arg,escript.Data):
296            return arg.getRank()
297        elif isinstance(arg,float):
298            return 0
299        elif isinstance(arg,int):
300            return 0
301        elif isinstance(arg,Symbol):
302            return arg.getRank()
303        else:
304          raise TypeError,"getShape: cannot identify shape"
305    def getShape(arg):
306      """      """
307      identifies the shape of its argument      identifies the shape of its argument
308    
# Line 291  def pokeShape(arg): Line 324  def pokeShape(arg):
324      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
325          return arg.getShape()          return arg.getShape()
326      else:      else:
327        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
328    
329  def pokeDim(arg):  def pokeDim(arg):
330      """      """
# Line 314  def commonShape(arg0,arg1): Line 347  def commonShape(arg0,arg1):
347      """      """
348      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.
349    
350      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
351      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
352      @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.
353      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
354      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
355      """      """
356      sh0=pokeShape(arg0)      sh0=getShape(arg0)
357      sh1=pokeShape(arg1)      sh1=getShape(arg1)
358      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
359         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
360               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 339  def commonDim(*args): Line 372  def commonDim(*args):
372      """      """
373      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.
374    
375      @param *args: given objects      @param args: given objects
376      @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
377               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
378      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 361  def testForZero(arg): Line 394  def testForZero(arg):
394    
395      @param arg: a given object      @param arg: a given object
396      @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}
397      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
398      @rtype : C{bool}      @rtype: C{bool}
399      """      """
400      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
401         return not Lsup(arg)>0.         return not Lsup(arg)>0.
# Line 395  def matchType(arg0=0.,arg1=0.): Line 428  def matchType(arg0=0.,arg1=0.):
428         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
429            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
430         elif isinstance(arg1,float):         elif isinstance(arg1,float):
431            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
432         elif isinstance(arg1,int):         elif isinstance(arg1,int):
433            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
434         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
435            pass            pass
436         else:         else:
# Line 421  def matchType(arg0=0.,arg1=0.): Line 454  def matchType(arg0=0.,arg1=0.):
454         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
455            pass            pass
456         elif isinstance(arg1,float):         elif isinstance(arg1,float):
457            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
458         elif isinstance(arg1,int):         elif isinstance(arg1,int):
459            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
460         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
461            pass            pass
462         else:         else:
463            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
464      elif isinstance(arg0,float):      elif isinstance(arg0,float):
465         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
466            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
467         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
468            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
469         elif isinstance(arg1,float):         elif isinstance(arg1,float):
470            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
471            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
472         elif isinstance(arg1,int):         elif isinstance(arg1,int):
473            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
474            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
475         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
476            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
477         else:         else:
478            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
479      elif isinstance(arg0,int):      elif isinstance(arg0,int):
480         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
481            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
482         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
483            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
484         elif isinstance(arg1,float):         elif isinstance(arg1,float):
485            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
486            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
487         elif isinstance(arg1,int):         elif isinstance(arg1,int):
488            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
489            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
490         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
491            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
492         else:         else:
493            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
494      else:      else:
# Line 465  def matchType(arg0=0.,arg1=0.): Line 498  def matchType(arg0=0.,arg1=0.):
498    
499  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
500      """      """
501            return representations of arg0 amd arg1 which ahve the same shape
   
     If shape is not given the shape "largest" shape of args is used.  
502    
503      @param args: a given ob      @param arg0: a given object
504      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
505      @return: True if the argument is identical to zero.      @param arg1: a given object
506      @rtype: C{list} of C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
507        @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
508        @rtype: C{tuple}
509      """      """
510      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
511      sh0=pokeShape(arg0)      sh0=getShape(arg0)
512      sh1=pokeShape(arg1)      sh1=getShape(arg1)
513      if len(sh0)<len(sh):      if len(sh0)<len(sh):
514         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
515      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
516         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
517      else:      else:
518         return arg0,arg1         return arg0,arg1
519  #=========================================================  #=========================================================
# Line 500  class Symbol(object): Line 533  class Symbol(object):
533         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
534         symbols or any other object.         symbols or any other object.
535    
536         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
537         @type arg: C{list}         @type args: C{list}
538         @param shape: the shape         @param shape: the shape
539         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
540         @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
# Line 544  class Symbol(object): Line 577  class Symbol(object):
577         """         """
578         the shape of the symbol.         the shape of the symbol.
579    
580         @return : the shape of the symbol.         @return: the shape of the symbol.
581         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
582         """         """
583         return self.__shape         return self.__shape
# Line 553  class Symbol(object): Line 586  class Symbol(object):
586         """         """
587         the spatial dimension         the spatial dimension
588    
589         @return : the spatial dimension         @return: the spatial dimension
590         @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.
591         """         """
592         return self.__dim         return self.__dim
# Line 577  class Symbol(object): Line 610  class Symbol(object):
610         """         """
611         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.
612    
613         @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].
614         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
615         @rtype: C{list} of objects         @rtype: C{list} of objects
616         @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.
617         """         """
618         out=[]         out=[]
619         for a in self.getArgument():         for a in self.getArgument():
# Line 604  class Symbol(object): Line 637  class Symbol(object):
637            if isinstance(a,Symbol):            if isinstance(a,Symbol):
638               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
639            else:            else:
640                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
641                if len(s)>0:                if len(s)>0:
642                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
643                else:                else:
644                   out.append(a)                   out.append(a)
645         return out         return out
# Line 696  class Symbol(object): Line 729  class Symbol(object):
729         else:         else:
730            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
731            if len(s)>0:            if len(s)>0:
732               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
733            else:            else:
734               return 0.               return 0.
735    
# Line 704  class Symbol(object): Line 737  class Symbol(object):
737         """         """
738         returns -self.         returns -self.
739    
740         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
741         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
742         """         """
743         return self*(-1.)         return self*(-1.)
# Line 713  class Symbol(object): Line 746  class Symbol(object):
746         """         """
747         returns +self.         returns +self.
748    
749         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
750         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
751         """         """
752         return self*(1.)         return self*(1.)
753    
754     def __abs__(self):     def __abs__(self):
755         """         """
756         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
757         """         """
758         return Abs_Symbol(self)         return Abs_Symbol(self)
759    
# Line 730  class Symbol(object): Line 763  class Symbol(object):
763    
764         @param other: object to be added to this object         @param other: object to be added to this object
765         @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}.
766         @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}
767         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
768         """         """
769         return add(self,other)         return add(self,other)
# Line 741  class Symbol(object): Line 774  class Symbol(object):
774    
775         @param other: object this object is added to         @param other: object this object is added to
776         @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}.
777         @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
778         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
779         """         """
780         return add(other,self)         return add(other,self)
# Line 752  class Symbol(object): Line 785  class Symbol(object):
785    
786         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
787         @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}.
788         @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
789         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
790         """         """
791         return add(self,-other)         return add(self,-other)
# Line 763  class Symbol(object): Line 796  class Symbol(object):
796    
797         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
798         @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}.
799         @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}.
800         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
801         """         """
802         return add(-self,other)         return add(-self,other)
# Line 774  class Symbol(object): Line 807  class Symbol(object):
807    
808         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
809         @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}.
810         @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}.
811         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
812         """         """
813         return mult(self,other)         return mult(self,other)
# Line 785  class Symbol(object): Line 818  class Symbol(object):
818    
819         @param other: object this object is multiplied with         @param other: object this object is multiplied with
820         @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}.
821         @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.
822         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
823         """         """
824         return mult(other,self)         return mult(other,self)
# Line 796  class Symbol(object): Line 829  class Symbol(object):
829    
830         @param other: object dividing this object         @param other: object dividing this object
831         @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}.
832         @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}
833         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
834         """         """
835         return quotient(self,other)         return quotient(self,other)
# Line 807  class Symbol(object): Line 840  class Symbol(object):
840    
841         @param other: object dividing this object         @param other: object dividing this object
842         @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}.
843         @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
844         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
845         """         """
846         return quotient(other,self)         return quotient(other,self)
# Line 818  class Symbol(object): Line 851  class Symbol(object):
851    
852         @param other: exponent         @param other: exponent
853         @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}.
854         @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}
855         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
856         """         """
857         return power(self,other)         return power(self,other)
# Line 829  class Symbol(object): Line 862  class Symbol(object):
862    
863         @param other: basis         @param other: basis
864         @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}.
865         @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
866         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
867         """         """
868         return power(other,self)         return power(other,self)
869    
870       def __getitem__(self,index):
871           """
872           returns the slice defined by index
873    
874           @param index: defines a
875           @type index: C{slice} or C{int} or a C{tuple} of them
876           @return: a L{Symbol} representing the slice defined by index
877           @rtype: L{DependendSymbol}
878           """
879           return GetSlice_Symbol(self,index)
880    
881  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
882     """     """
883     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.
884     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  
885        
886     Example:     Example::
887        
888     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
889     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
890     print u1==u2       print u1==u2
891     False       False
892        
893        but       but::
894    
895     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
896     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
897     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
898     print u1==u2, u1==u3       print u1==u2, u1==u3
899     True False       True False
900    
901     @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.
902     """     """
# Line 884  class DependendSymbol(Symbol): Line 928  class DependendSymbol(Symbol):
928  #=========================================================  #=========================================================
929  #  Unary operations prserving the shape  #  Unary operations prserving the shape
930  #========================================================  #========================================================
931    class GetSlice_Symbol(DependendSymbol):
932       """
933       L{Symbol} representing getting a slice for a L{Symbol}
934       """
935       def __init__(self,arg,index):
936          """
937          initialization of wherePositive L{Symbol} with argument arg
938          @param arg: argument
939          @type arg: L{Symbol}.
940          @param index: defines index
941          @type index: C{slice} or C{int} or a C{tuple} of them
942          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
943          @raises ValueError: if a step is given
944          """
945          if not isinstance(index,tuple): index=(index,)
946          if len(index)>arg.getRank():
947               raise IndexError,"GetSlice_Symbol: index out of range."
948          sh=()
949          index2=()
950          for i in range(len(index)):
951             ix=index[i]
952             if isinstance(ix,int):
953                if ix<0 or ix>=arg.getShape()[i]:
954                   raise ValueError,"GetSlice_Symbol: index out of range."
955                index2=index2+(ix,)
956             else:
957               if not ix.step==None:
958                 raise ValueError,"GetSlice_Symbol: steping is not supported."
959               if ix.start==None:
960                  s=0
961               else:
962                  s=ix.start
963               if ix.stop==None:
964                  e=arg.getShape()[i]
965               else:
966                  e=ix.stop
967                  if e>arg.getShape()[i]:
968                     raise IndexError,"GetSlice_Symbol: index out of range."
969               index2=index2+(slice(s,e),)
970               if e>s:
971                   sh=sh+(e-s,)
972               elif s>e:
973                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
974          for i in range(len(index),arg.getRank()):
975              index2=index2+(slice(0,arg.getShape()[i]),)
976              sh=sh+(arg.getShape()[i],)
977          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
978    
979       def getMyCode(self,argstrs,format="escript"):
980          """
981          returns a program code that can be used to evaluate the symbol.
982    
983          @param argstrs: gives for each argument a string representing the argument for the evaluation.
984          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
985          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
986          @type format: C{str}
987          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
988          @rtype: C{str}
989          @raise NotImplementedError: if the requested format is not available
990          """
991          if format=="escript" or format=="str"  or format=="text":
992             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
993          else:
994             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
995    
996       def substitute(self,argvals):
997          """
998          assigns new values to symbols in the definition of the symbol.
999          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1000    
1001          @param argvals: new values assigned to symbols
1002          @type argvals: C{dict} with keywords of type L{Symbol}.
1003          @return: result of the substitution process. Operations are executed as much as possible.
1004          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1005          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1006          """
1007          if argvals.has_key(self):
1008             arg=argvals[self]
1009             if self.isAppropriateValue(arg):
1010                return arg
1011             else:
1012                raise TypeError,"%s: new value is not appropriate."%str(self)
1013          else:
1014             args=self.getSubstitutedArguments(argvals)
1015             arg=args[0]
1016             index=args[1]
1017             return arg.__getitem__(index)
1018    
1019  def log10(arg):  def log10(arg):
1020     """     """
1021     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
1022    
1023     @param arg: argument     @param arg: argument
1024     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1025     @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.
1026     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1027     """     """
1028     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 912  def wherePositive(arg): Line 1044  def wherePositive(arg):
1044    
1045     @param arg: argument     @param arg: argument
1046     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1047     @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.
1048     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1049     """     """
1050     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1051        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1052        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1053        return out        return out
1054     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1055        return arg._wherePositive()        return arg._wherePositive()
# Line 958  class WherePositive_Symbol(DependendSymb Line 1090  class WherePositive_Symbol(DependendSymb
1090        @type format: C{str}        @type format: C{str}
1091        @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.
1092        @rtype: C{str}        @rtype: C{str}
1093        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1094        """        """
1095        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1096            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 994  def whereNegative(arg): Line 1126  def whereNegative(arg):
1126    
1127     @param arg: argument     @param arg: argument
1128     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1129     @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.
1130     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1131     """     """
1132     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1133        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1134        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1135        return out        return out
1136     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1137        return arg._whereNegative()        return arg._whereNegative()
# Line 1040  class WhereNegative_Symbol(DependendSymb Line 1172  class WhereNegative_Symbol(DependendSymb
1172        @type format: C{str}        @type format: C{str}
1173        @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.
1174        @rtype: C{str}        @rtype: C{str}
1175        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1176        """        """
1177        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1178            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1076  def whereNonNegative(arg): Line 1208  def whereNonNegative(arg):
1208    
1209     @param arg: argument     @param arg: argument
1210     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1211     @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.
1212     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1213     """     """
1214     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1215        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1216        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1217        return out        return out
1218     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1219        return arg._whereNonNegative()        return arg._whereNonNegative()
# Line 1106  def whereNonPositive(arg): Line 1238  def whereNonPositive(arg):
1238    
1239     @param arg: argument     @param arg: argument
1240     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1241     @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.
1242     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1243     """     """
1244     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1245        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1246        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1247        return out        return out
1248     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1249        return arg._whereNonPositive()        return arg._whereNonPositive()
# Line 1138  def whereZero(arg,tol=0.): Line 1270  def whereZero(arg,tol=0.):
1270     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1271     @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.
1272     @type tol: C{float}     @type tol: C{float}
1273     @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.
1274     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1275     """     """
1276     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1277        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1278        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1279        return out        return out
1280     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1281        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1282     elif isinstance(arg,float):     elif isinstance(arg,float):
1283        if abs(arg)<=tol:        if abs(arg)<=tol:
1284          return 1.          return 1.
# Line 1187  class WhereZero_Symbol(DependendSymbol): Line 1316  class WhereZero_Symbol(DependendSymbol):
1316        @type format: C{str}        @type format: C{str}
1317        @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.
1318        @rtype: C{str}        @rtype: C{str}
1319        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1320        """        """
1321        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1322           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1221  def whereNonZero(arg,tol=0.): Line 1350  def whereNonZero(arg,tol=0.):
1350    
1351     @param arg: argument     @param arg: argument
1352     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1353     @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.
1354     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1355     """     """
1356     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1357        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1358        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1359        return out        return out
1360     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1361        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1362     elif isinstance(arg,float):     elif isinstance(arg,float):
1363        if abs(arg)>tol:        if abs(arg)>tol:
1364          return 1.          return 1.
# Line 1248  def whereNonZero(arg,tol=0.): Line 1374  def whereNonZero(arg,tol=0.):
1374     else:     else:
1375        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1376    
1377    def erf(arg):
1378       """
1379       returns erf of argument arg
1380    
1381       @param arg: argument
1382       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1383       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1384       @raises TypeError: if the type of the argument is not expected.
1385       """
1386       if isinstance(arg,escript.Data):
1387          return arg._erf()
1388       else:
1389          raise TypeError,"erf: Unknown argument type."
1390    
1391  def sin(arg):  def sin(arg):
1392     """     """
1393     returns sine of argument arg     returns sine of argument arg
1394    
1395     @param arg: argument     @param arg: argument
1396     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1397     @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.
1398     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1399     """     """
1400     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1292  class Sin_Symbol(DependendSymbol): Line 1432  class Sin_Symbol(DependendSymbol):
1432        @type format: C{str}        @type format: C{str}
1433        @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.
1434        @rtype: C{str}        @rtype: C{str}
1435        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1436        """        """
1437        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1438            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1344  def cos(arg): Line 1484  def cos(arg):
1484    
1485     @param arg: argument     @param arg: argument
1486     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1487     @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.
1488     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1489     """     """
1490     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1382  class Cos_Symbol(DependendSymbol): Line 1522  class Cos_Symbol(DependendSymbol):
1522        @type format: C{str}        @type format: C{str}
1523        @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.
1524        @rtype: C{str}        @rtype: C{str}
1525        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1526        """        """
1527        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1528            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1434  def tan(arg): Line 1574  def tan(arg):
1574    
1575     @param arg: argument     @param arg: argument
1576     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1577     @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.
1578     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1579     """     """
1580     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1472  class Tan_Symbol(DependendSymbol): Line 1612  class Tan_Symbol(DependendSymbol):
1612        @type format: C{str}        @type format: C{str}
1613        @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.
1614        @rtype: C{str}        @rtype: C{str}
1615        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1616        """        """
1617        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1618            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1524  def asin(arg): Line 1664  def asin(arg):
1664    
1665     @param arg: argument     @param arg: argument
1666     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1667     @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.
1668     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1669     """     """
1670     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1562  class Asin_Symbol(DependendSymbol): Line 1702  class Asin_Symbol(DependendSymbol):
1702        @type format: C{str}        @type format: C{str}
1703        @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.
1704        @rtype: C{str}        @rtype: C{str}
1705        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1706        """        """
1707        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1708            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1614  def acos(arg): Line 1754  def acos(arg):
1754    
1755     @param arg: argument     @param arg: argument
1756     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1757     @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.
1758     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1759     """     """
1760     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1652  class Acos_Symbol(DependendSymbol): Line 1792  class Acos_Symbol(DependendSymbol):
1792        @type format: C{str}        @type format: C{str}
1793        @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.
1794        @rtype: C{str}        @rtype: C{str}
1795        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1796        """        """
1797        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1798            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1704  def atan(arg): Line 1844  def atan(arg):
1844    
1845     @param arg: argument     @param arg: argument
1846     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1847     @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.
1848     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1849     """     """
1850     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1742  class Atan_Symbol(DependendSymbol): Line 1882  class Atan_Symbol(DependendSymbol):
1882        @type format: C{str}        @type format: C{str}
1883        @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.
1884        @rtype: C{str}        @rtype: C{str}
1885        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1886        """        """
1887        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1888            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1794  def sinh(arg): Line 1934  def sinh(arg):
1934    
1935     @param arg: argument     @param arg: argument
1936     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1937     @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.
1938     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1939     """     """
1940     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1832  class Sinh_Symbol(DependendSymbol): Line 1972  class Sinh_Symbol(DependendSymbol):
1972        @type format: C{str}        @type format: C{str}
1973        @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.
1974        @rtype: C{str}        @rtype: C{str}
1975        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1976        """        """
1977        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1978            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1884  def cosh(arg): Line 2024  def cosh(arg):
2024    
2025     @param arg: argument     @param arg: argument
2026     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2027     @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.
2028     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2029     """     """
2030     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1922  class Cosh_Symbol(DependendSymbol): Line 2062  class Cosh_Symbol(DependendSymbol):
2062        @type format: C{str}        @type format: C{str}
2063        @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.
2064        @rtype: C{str}        @rtype: C{str}
2065        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2066        """        """
2067        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2068            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1974  def tanh(arg): Line 2114  def tanh(arg):
2114    
2115     @param arg: argument     @param arg: argument
2116     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2117     @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.
2118     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2119     """     """
2120     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2012  class Tanh_Symbol(DependendSymbol): Line 2152  class Tanh_Symbol(DependendSymbol):
2152        @type format: C{str}        @type format: C{str}
2153        @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.
2154        @rtype: C{str}        @rtype: C{str}
2155        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2156        """        """
2157        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2158            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2064  def asinh(arg): Line 2204  def asinh(arg):
2204    
2205     @param arg: argument     @param arg: argument
2206     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2207     @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.
2208     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2209     """     """
2210     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2102  class Asinh_Symbol(DependendSymbol): Line 2242  class Asinh_Symbol(DependendSymbol):
2242        @type format: C{str}        @type format: C{str}
2243        @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.
2244        @rtype: C{str}        @rtype: C{str}
2245        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2246        """        """
2247        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2248            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2154  def acosh(arg): Line 2294  def acosh(arg):
2294    
2295     @param arg: argument     @param arg: argument
2296     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2297     @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.
2298     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2299     """     """
2300     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2192  class Acosh_Symbol(DependendSymbol): Line 2332  class Acosh_Symbol(DependendSymbol):
2332        @type format: C{str}        @type format: C{str}
2333        @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.
2334        @rtype: C{str}        @rtype: C{str}
2335        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2336        """        """
2337        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2338            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2244  def atanh(arg): Line 2384  def atanh(arg):
2384    
2385     @param arg: argument     @param arg: argument
2386     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2387     @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.
2388     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2389     """     """
2390     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2282  class Atanh_Symbol(DependendSymbol): Line 2422  class Atanh_Symbol(DependendSymbol):
2422        @type format: C{str}        @type format: C{str}
2423        @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.
2424        @rtype: C{str}        @rtype: C{str}
2425        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2426        """        """
2427        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2428            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2334  def exp(arg): Line 2474  def exp(arg):
2474    
2475     @param arg: argument     @param arg: argument
2476     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2477     @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.
2478     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2479     """     """
2480     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2372  class Exp_Symbol(DependendSymbol): Line 2512  class Exp_Symbol(DependendSymbol):
2512        @type format: C{str}        @type format: C{str}
2513        @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.
2514        @rtype: C{str}        @rtype: C{str}
2515        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2516        """        """
2517        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2518            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2424  def sqrt(arg): Line 2564  def sqrt(arg):
2564    
2565     @param arg: argument     @param arg: argument
2566     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2567     @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.
2568     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2569     """     """
2570     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2462  class Sqrt_Symbol(DependendSymbol): Line 2602  class Sqrt_Symbol(DependendSymbol):
2602        @type format: C{str}        @type format: C{str}
2603        @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.
2604        @rtype: C{str}        @rtype: C{str}
2605        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2606        """        """
2607        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2608            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2514  def log(arg): Line 2654  def log(arg):
2654    
2655     @param arg: argument     @param arg: argument
2656     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2657     @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.
2658     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2659     """     """
2660     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2552  class Log_Symbol(DependendSymbol): Line 2692  class Log_Symbol(DependendSymbol):
2692        @type format: C{str}        @type format: C{str}
2693        @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.
2694        @rtype: C{str}        @rtype: C{str}
2695        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2696        """        """
2697        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2698            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2604  def sign(arg): Line 2744  def sign(arg):
2744    
2745     @param arg: argument     @param arg: argument
2746     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2747     @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.
2748     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2749     """     """
2750     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2652  class Abs_Symbol(DependendSymbol): Line 2792  class Abs_Symbol(DependendSymbol):
2792        @type format: C{str}        @type format: C{str}
2793        @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.
2794        @rtype: C{str}        @rtype: C{str}
2795        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2796        """        """
2797        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2798            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2704  def minval(arg): Line 2844  def minval(arg):
2844    
2845     @param arg: argument     @param arg: argument
2846     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2847     @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.
2848     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2849     """     """
2850     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2745  class Minval_Symbol(DependendSymbol): Line 2885  class Minval_Symbol(DependendSymbol):
2885        @type format: C{str}        @type format: C{str}
2886        @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.
2887        @rtype: C{str}        @rtype: C{str}
2888        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2889        """        """
2890        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2891            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2781  def maxval(arg): Line 2921  def maxval(arg):
2921    
2922     @param arg: argument     @param arg: argument
2923     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2924     @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.
2925     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2926     """     """
2927     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2822  class Maxval_Symbol(DependendSymbol): Line 2962  class Maxval_Symbol(DependendSymbol):
2962        @type format: C{str}        @type format: C{str}
2963        @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.
2964        @rtype: C{str}        @rtype: C{str}
2965        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2966        """        """
2967        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2968            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2858  def length(arg): Line 2998  def length(arg):
2998    
2999     @param arg: argument     @param arg: argument
3000     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3001     @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.
3002     """     """
3003     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
3004    
3005    def trace(arg,axis_offset=0):
3006       """
3007       returns the trace of arg which the sum of arg[k,k] over k.
3008    
3009       @param arg: argument
3010       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3011       @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
3012                      C{axis_offset} and axis_offset+1 must be equal.
3013       @type axis_offset: C{int}
3014       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
3015       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3016       """
3017       if isinstance(arg,numarray.NumArray):
3018          sh=arg.shape
3019          if len(sh)<2:
3020            raise ValueError,"rank of argument must be greater than 1"
3021          if axis_offset<0 or axis_offset>len(sh)-2:
3022            raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
3023          s1=1
3024          for i in range(axis_offset): s1*=sh[i]
3025          s2=1
3026          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3027          if not sh[axis_offset] == sh[axis_offset+1]:
3028            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3029          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
3030          out=numarray.zeros([s1,s2],numarray.Float64)
3031          for i1 in range(s1):
3032            for i2 in range(s2):
3033                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
3034          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3035          return out
3036       elif isinstance(arg,escript.Data):
3037          if arg.getRank()<2:
3038            raise ValueError,"rank of argument must be greater than 1"
3039          if axis_offset<0 or axis_offset>arg.getRank()-2:
3040            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3041          s=list(arg.getShape())        
3042          if not s[axis_offset] == s[axis_offset+1]:
3043            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3044          return arg._trace(axis_offset)
3045       elif isinstance(arg,float):
3046          raise TypeError,"illegal argument type float."
3047       elif isinstance(arg,int):
3048          raise TypeError,"illegal argument type int."
3049       elif isinstance(arg,Symbol):
3050          return Trace_Symbol(arg,axis_offset)
3051       else:
3052          raise TypeError,"Unknown argument type."
3053    
3054    class Trace_Symbol(DependendSymbol):
3055       """
3056       L{Symbol} representing the result of the trace function
3057       """
3058       def __init__(self,arg,axis_offset=0):
3059          """
3060          initialization of trace L{Symbol} with argument arg
3061          @param arg: argument of function
3062          @type arg: L{Symbol}.
3063          @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
3064                      C{axis_offset} and axis_offset+1 must be equal.
3065          @type axis_offset: C{int}
3066          """
3067          if arg.getRank()<2:
3068            raise ValueError,"rank of argument must be greater than 1"
3069          if axis_offset<0 or axis_offset>arg.getRank()-2:
3070            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3071          s=list(arg.getShape())        
3072          if not s[axis_offset] == s[axis_offset+1]:
3073            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3074          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3075    
3076       def getMyCode(self,argstrs,format="escript"):
3077          """
3078          returns a program code that can be used to evaluate the symbol.
3079    
3080          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3081          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3082          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3083          @type format: C{str}
3084          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3085          @rtype: C{str}
3086          @raise NotImplementedError: if the requested format is not available
3087          """
3088          if format=="escript" or format=="str"  or format=="text":
3089             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3090          else:
3091             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3092    
3093       def substitute(self,argvals):
3094          """
3095          assigns new values to symbols in the definition of the symbol.
3096          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3097    
3098          @param argvals: new values assigned to symbols
3099          @type argvals: C{dict} with keywords of type L{Symbol}.
3100          @return: result of the substitution process. Operations are executed as much as possible.
3101          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3102          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3103          """
3104          if argvals.has_key(self):
3105             arg=argvals[self]
3106             if self.isAppropriateValue(arg):
3107                return arg
3108             else:
3109                raise TypeError,"%s: new value is not appropriate."%str(self)
3110          else:
3111             arg=self.getSubstitutedArguments(argvals)
3112             return trace(arg[0],axis_offset=arg[1])
3113    
3114       def diff(self,arg):
3115          """
3116          differential of this object
3117    
3118          @param arg: the derivative is calculated with respect to arg
3119          @type arg: L{escript.Symbol}
3120          @return: derivative with respect to C{arg}
3121          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3122          """
3123          if arg==self:
3124             return identity(self.getShape())
3125          else:
3126             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3127    
3128    def transpose(arg,axis_offset=None):
3129       """
3130       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3131    
3132       @param arg: argument
3133       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3134       @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.
3135                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3136       @type axis_offset: C{int}
3137       @return: transpose of arg
3138       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3139       """
3140       if isinstance(arg,numarray.NumArray):
3141          if axis_offset==None: axis_offset=int(arg.rank/2)
3142          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3143       elif isinstance(arg,escript.Data):
3144          r=arg.getRank()
3145          if axis_offset==None: axis_offset=int(r/2)
3146          if axis_offset<0 or axis_offset>r:
3147            raise ValueError,"axis_offset must be between 0 and %s"%r
3148          return arg._transpose(axis_offset)
3149       elif isinstance(arg,float):
3150          if not ( axis_offset==0 or axis_offset==None):
3151            raise ValueError,"axis_offset must be 0 for float argument"
3152          return arg
3153       elif isinstance(arg,int):
3154          if not ( axis_offset==0 or axis_offset==None):
3155            raise ValueError,"axis_offset must be 0 for int argument"
3156          return float(arg)
3157       elif isinstance(arg,Symbol):
3158          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3159          return Transpose_Symbol(arg,axis_offset)
3160       else:
3161          raise TypeError,"Unknown argument type."
3162    
3163    class Transpose_Symbol(DependendSymbol):
3164       """
3165       L{Symbol} representing the result of the transpose function
3166       """
3167       def __init__(self,arg,axis_offset=None):
3168          """
3169          initialization of transpose L{Symbol} with argument arg
3170    
3171          @param arg: argument of function
3172          @type arg: L{Symbol}.
3173          @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.
3174                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3175          @type axis_offset: C{int}
3176          """
3177          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3178          if axis_offset<0 or axis_offset>arg.getRank():
3179            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3180          s=arg.getShape()
3181          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3182    
3183       def getMyCode(self,argstrs,format="escript"):
3184          """
3185          returns a program code that can be used to evaluate the symbol.
3186    
3187          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3188          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3189          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3190          @type format: C{str}
3191          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3192          @rtype: C{str}
3193          @raise NotImplementedError: if the requested format is not available
3194          """
3195          if format=="escript" or format=="str"  or format=="text":
3196             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3197          else:
3198             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3199    
3200       def substitute(self,argvals):
3201          """
3202          assigns new values to symbols in the definition of the symbol.
3203          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3204    
3205          @param argvals: new values assigned to symbols
3206          @type argvals: C{dict} with keywords of type L{Symbol}.
3207          @return: result of the substitution process. Operations are executed as much as possible.
3208          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3209          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3210          """
3211          if argvals.has_key(self):
3212             arg=argvals[self]
3213             if self.isAppropriateValue(arg):
3214                return arg
3215             else:
3216                raise TypeError,"%s: new value is not appropriate."%str(self)
3217          else:
3218             arg=self.getSubstitutedArguments(argvals)
3219             return transpose(arg[0],axis_offset=arg[1])
3220    
3221       def diff(self,arg):
3222          """
3223          differential of this object
3224    
3225          @param arg: the derivative is calculated with respect to arg
3226          @type arg: L{escript.Symbol}
3227          @return: derivative with respect to C{arg}
3228          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3229          """
3230          if arg==self:
3231             return identity(self.getShape())
3232          else:
3233             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3234    
3235    def swap_axes(arg,axis0=0,axis1=1):
3236       """
3237       returns the swap of arg by swaping the components axis0 and axis1
3238    
3239       @param arg: argument
3240       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3241       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3242       @type axis0: C{int}
3243       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3244       @type axis1: C{int}
3245       @return: C{arg} with swaped components
3246       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3247       """
3248       if axis0 > axis1:
3249          axis0,axis1=axis1,axis0
3250       if isinstance(arg,numarray.NumArray):
3251          return numarray.swapaxes(arg,axis0,axis1)
3252       elif isinstance(arg,escript.Data):
3253          return arg._swap_axes(axis0,axis1)
3254       elif isinstance(arg,float):
3255          raise TyepError,"float argument is not supported."
3256       elif isinstance(arg,int):
3257          raise TyepError,"int argument is not supported."
3258       elif isinstance(arg,Symbol):
3259          return SwapAxes_Symbol(arg,axis0,axis1)
3260       else:
3261          raise TypeError,"Unknown argument type."
3262    
3263    class SwapAxes_Symbol(DependendSymbol):
3264       """
3265       L{Symbol} representing the result of the swap function
3266       """
3267       def __init__(self,arg,axis0=0,axis1=1):
3268          """
3269          initialization of swap L{Symbol} with argument arg
3270    
3271          @param arg: argument
3272          @type arg: L{Symbol}.
3273          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3274          @type axis0: C{int}
3275          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3276          @type axis1: C{int}
3277          """
3278          if arg.getRank()<2:
3279             raise ValueError,"argument must have at least rank 2."
3280          if axis0<0 or axis0>arg.getRank()-1:
3281             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3282          if axis1<0 or axis1>arg.getRank()-1:
3283             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3284          if axis0 == axis1:
3285             raise ValueError,"axis indices must be different."
3286          if axis0 > axis1:
3287             axis0,axis1=axis1,axis0
3288          s=arg.getShape()
3289          s_out=[]
3290          for i in range(len(s)):
3291             if i == axis0:
3292                s_out.append(s[axis1])
3293             elif i == axis1:
3294                s_out.append(s[axis0])
3295             else:
3296                s_out.append(s[i])
3297          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3298    
3299       def getMyCode(self,argstrs,format="escript"):
3300          """
3301          returns a program code that can be used to evaluate the symbol.
3302    
3303          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3304          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3305          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3306          @type format: C{str}
3307          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3308          @rtype: C{str}
3309          @raise NotImplementedError: if the requested format is not available
3310          """
3311          if format=="escript" or format=="str"  or format=="text":
3312             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3313          else:
3314             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3315    
3316       def substitute(self,argvals):
3317          """
3318          assigns new values to symbols in the definition of the symbol.
3319          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3320    
3321          @param argvals: new values assigned to symbols
3322          @type argvals: C{dict} with keywords of type L{Symbol}.
3323          @return: result of the substitution process. Operations are executed as much as possible.
3324          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3325          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3326          """
3327          if argvals.has_key(self):
3328             arg=argvals[self]
3329             if self.isAppropriateValue(arg):
3330                return arg
3331             else:
3332                raise TypeError,"%s: new value is not appropriate."%str(self)
3333          else:
3334             arg=self.getSubstitutedArguments(argvals)
3335             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3336    
3337       def diff(self,arg):
3338          """
3339          differential of this object
3340    
3341          @param arg: the derivative is calculated with respect to arg
3342          @type arg: L{escript.Symbol}
3343          @return: derivative with respect to C{arg}
3344          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3345          """
3346          if arg==self:
3347             return identity(self.getShape())
3348          else:
3349             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3350    
3351    def symmetric(arg):
3352        """
3353        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3354    
3355        @param arg: square matrix. Must have rank 2 or 4 and be square.
3356        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3357        @return: symmetric part of arg
3358        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3359        """
3360        if isinstance(arg,numarray.NumArray):
3361          if arg.rank==2:
3362            if not (arg.shape[0]==arg.shape[1]):
3363               raise ValueError,"argument must be square."
3364          elif arg.rank==4:
3365            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3366               raise ValueError,"argument must be square."
3367          else:
3368            raise ValueError,"rank 2 or 4 is required."
3369          return (arg+transpose(arg))/2
3370        elif isinstance(arg,escript.Data):
3371          if arg.getRank()==2:
3372            if not (arg.getShape()[0]==arg.getShape()[1]):
3373               raise ValueError,"argument must be square."
3374            return arg._symmetric()
3375          elif arg.getRank()==4:
3376            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3377               raise ValueError,"argument must be square."
3378            return arg._symmetric()
3379          else:
3380            raise ValueError,"rank 2 or 4 is required."
3381        elif isinstance(arg,float):
3382          return arg
3383        elif isinstance(arg,int):
3384          return float(arg)
3385        elif isinstance(arg,Symbol):
3386          if arg.getRank()==2:
3387            if not (arg.getShape()[0]==arg.getShape()[1]):
3388               raise ValueError,"argument must be square."
3389          elif arg.getRank()==4:
3390            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3391               raise ValueError,"argument must be square."
3392          else:
3393            raise ValueError,"rank 2 or 4 is required."
3394          return (arg+transpose(arg))/2
3395        else:
3396          raise TypeError,"symmetric: Unknown argument type."
3397    
3398    def nonsymmetric(arg):
3399        """
3400        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3401    
3402        @param arg: square matrix. Must have rank 2 or 4 and be square.
3403        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3404        @return: nonsymmetric part of arg
3405        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3406        """
3407        if isinstance(arg,numarray.NumArray):
3408          if arg.rank==2:
3409            if not (arg.shape[0]==arg.shape[1]):
3410               raise ValueError,"nonsymmetric: argument must be square."
3411          elif arg.rank==4:
3412            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3413               raise ValueError,"nonsymmetric: argument must be square."
3414          else:
3415            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3416          return (arg-transpose(arg))/2
3417        elif isinstance(arg,escript.Data):
3418          if arg.getRank()==2:
3419            if not (arg.getShape()[0]==arg.getShape()[1]):
3420               raise ValueError,"argument must be square."
3421            return arg._nonsymmetric()
3422          elif arg.getRank()==4:
3423            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3424               raise ValueError,"argument must be square."
3425            return arg._nonsymmetric()
3426          else:
3427            raise ValueError,"rank 2 or 4 is required."
3428        elif isinstance(arg,float):
3429          return arg
3430        elif isinstance(arg,int):
3431          return float(arg)
3432        elif isinstance(arg,Symbol):
3433          if arg.getRank()==2:
3434            if not (arg.getShape()[0]==arg.getShape()[1]):
3435               raise ValueError,"nonsymmetric: argument must be square."
3436          elif arg.getRank()==4:
3437            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3438               raise ValueError,"nonsymmetric: argument must be square."
3439          else:
3440            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3441          return (arg-transpose(arg))/2
3442        else:
3443          raise TypeError,"nonsymmetric: Unknown argument type."
3444    
3445    def inverse(arg):
3446        """
3447        returns the inverse of the square matrix arg.
3448    
3449        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3450        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3451        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3452        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3453        @note: for L{escript.Data} objects the dimension is restricted to 3.
3454        """
3455        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3456        if isinstance(arg,numarray.NumArray):
3457          return numarray.linear_algebra.inverse(arg)
3458        elif isinstance(arg,escript.Data):
3459          return escript_inverse(arg)
3460        elif isinstance(arg,float):
3461          return 1./arg
3462        elif isinstance(arg,int):
3463          return 1./float(arg)
3464        elif isinstance(arg,Symbol):
3465          return Inverse_Symbol(arg)
3466        else:
3467          raise TypeError,"inverse: Unknown argument type."
3468    
3469    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3470          "arg is a Data objects!!!"
3471          if not arg.getRank()==2:
3472            raise ValueError,"escript_inverse: argument must have rank 2"
3473          s=arg.getShape()      
3474          if not s[0] == s[1]:
3475            raise ValueError,"escript_inverse: argument must be a square matrix."
3476          out=escript.Data(0.,s,arg.getFunctionSpace())
3477          if s[0]==1:
3478              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3479                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3480              out[0,0]=1./arg[0,0]
3481          elif s[0]==2:
3482              A11=arg[0,0]
3483              A12=arg[0,1]
3484              A21=arg[1,0]
3485              A22=arg[1,1]
3486              D = A11*A22-A12*A21
3487              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3488                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3489              D=1./D
3490              out[0,0]= A22*D
3491              out[1,0]=-A21*D
3492              out[0,1]=-A12*D
3493              out[1,1]= A11*D
3494          elif s[0]==3:
3495              A11=arg[0,0]
3496              A21=arg[1,0]
3497              A31=arg[2,0]
3498              A12=arg[0,1]
3499              A22=arg[1,1]
3500              A32=arg[2,1]
3501              A13=arg[0,2]
3502              A23=arg[1,2]
3503              A33=arg[2,2]
3504              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3505              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3506                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3507              D=1./D
3508              out[0,0]=(A22*A33-A23*A32)*D
3509              out[1,0]=(A31*A23-A21*A33)*D
3510              out[2,0]=(A21*A32-A31*A22)*D
3511              out[0,1]=(A13*A32-A12*A33)*D
3512              out[1,1]=(A11*A33-A31*A13)*D
3513              out[2,1]=(A12*A31-A11*A32)*D
3514              out[0,2]=(A12*A23-A13*A22)*D
3515              out[1,2]=(A13*A21-A11*A23)*D
3516              out[2,2]=(A11*A22-A12*A21)*D
3517          else:
3518             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3519          return out
3520    
3521    class Inverse_Symbol(DependendSymbol):
3522       """
3523       L{Symbol} representing the result of the inverse function
3524       """
3525       def __init__(self,arg):
3526          """
3527          initialization of inverse L{Symbol} with argument arg
3528          @param arg: argument of function
3529          @type arg: L{Symbol}.
3530          """
3531          if not arg.getRank()==2:
3532            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3533          s=arg.getShape()
3534          if not s[0] == s[1]:
3535            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3536          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3537    
3538       def getMyCode(self,argstrs,format="escript"):
3539          """
3540          returns a program code that can be used to evaluate the symbol.
3541    
3542          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3543          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3544          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3545          @type format: C{str}
3546          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3547          @rtype: C{str}
3548          @raise NotImplementedError: if the requested format is not available
3549          """
3550          if format=="escript" or format=="str"  or format=="text":
3551             return "inverse(%s)"%argstrs[0]
3552          else:
3553             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3554    
3555       def substitute(self,argvals):
3556          """
3557          assigns new values to symbols in the definition of the symbol.
3558          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3559    
3560          @param argvals: new values assigned to symbols
3561          @type argvals: C{dict} with keywords of type L{Symbol}.
3562          @return: result of the substitution process. Operations are executed as much as possible.
3563          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3564          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3565          """
3566          if argvals.has_key(self):
3567             arg=argvals[self]
3568             if self.isAppropriateValue(arg):
3569                return arg
3570             else:
3571                raise TypeError,"%s: new value is not appropriate."%str(self)
3572          else:
3573             arg=self.getSubstitutedArguments(argvals)
3574             return inverse(arg[0])
3575    
3576       def diff(self,arg):
3577          """
3578          differential of this object
3579    
3580          @param arg: the derivative is calculated with respect to arg
3581          @type arg: L{escript.Symbol}
3582          @return: derivative with respect to C{arg}
3583          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3584          """
3585          if arg==self:
3586             return identity(self.getShape())
3587          else:
3588             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3589    
3590    def eigenvalues(arg):
3591        """
3592        returns the eigenvalues of the square matrix arg.
3593    
3594        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3595                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3596        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3597        @return: the eigenvalues in increasing order.
3598        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3599        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3600        """
3601        if isinstance(arg,numarray.NumArray):
3602          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3603          out.sort()
3604          return out
3605        elif isinstance(arg,escript.Data):
3606          return arg._eigenvalues()
3607        elif isinstance(arg,Symbol):
3608          if not arg.getRank()==2:
3609            raise ValueError,"eigenvalues: argument must have rank 2"
3610          s=arg.getShape()      
3611          if not s[0] == s[1]:
3612            raise ValueError,"eigenvalues: argument must be a square matrix."
3613          if s[0]==1:
3614              return arg[0]
3615          elif s[0]==2:
3616              arg1=symmetric(arg)
3617              A11=arg1[0,0]
3618              A12=arg1[0,1]
3619              A22=arg1[1,1]
3620              trA=(A11+A22)/2.
3621              A11-=trA
3622              A22-=trA
3623              s=sqrt(A12**2-A11*A22)
3624              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3625          elif s[0]==3:
3626              arg1=symmetric(arg)
3627              A11=arg1[0,0]
3628              A12=arg1[0,1]
3629              A22=arg1[1,1]
3630              A13=arg1[0,2]
3631              A23=arg1[1,2]
3632              A33=arg1[2,2]
3633              trA=(A11+A22+A33)/3.
3634              A11-=trA
3635              A22-=trA
3636              A33-=trA
3637              A13_2=A13**2
3638              A23_2=A23**2
3639              A12_2=A12**2
3640              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3641              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3642              sq_p=sqrt(p/3.)
3643              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
3644              sq_p*=2.
3645              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3646               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3647               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3648              return trA+sq_p*f
3649          else:
3650             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3651        elif isinstance(arg,float):
3652          return arg
3653        elif isinstance(arg,int):
3654          return float(arg)
3655        else:
3656          raise TypeError,"eigenvalues: Unknown argument type."
3657    
3658    def eigenvalues_and_eigenvectors(arg):
3659        """
3660        returns the eigenvalues and eigenvectors of the square matrix arg.
3661    
3662        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3663                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3664        @type arg: L{escript.Data}
3665        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3666                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3667                 the eigenvector coresponding to the i-th eigenvalue.
3668        @rtype: L{tuple} of L{escript.Data}.
3669        @note: The dimension is restricted to 3.
3670        """
3671        if isinstance(arg,numarray.NumArray):
3672          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3673        elif isinstance(arg,escript.Data):
3674          return arg._eigenvalues_and_eigenvectors()
3675        elif isinstance(arg,Symbol):
3676          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3677        elif isinstance(arg,float):
3678          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3679        elif isinstance(arg,int):
3680          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3681        else:
3682          raise TypeError,"eigenvalues: Unknown argument type."
3683  #=======================================================  #=======================================================
3684  #  Binary operations:  #  Binary operations:
3685  #=======================================================  #=======================================================
# Line 2905  class Add_Symbol(DependendSymbol): Line 3723  class Add_Symbol(DependendSymbol):
3723         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3724         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3725         """         """
3726         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3727         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3728         if not sh0==sh1:         if not sh0==sh1:
3729            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3730         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 2921  class Add_Symbol(DependendSymbol): Line 3739  class Add_Symbol(DependendSymbol):
3739        @type format: C{str}        @type format: C{str}
3740        @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.
3741        @rtype: C{str}        @rtype: C{str}
3742        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3743        """        """
3744        if format=="str" or format=="text":        if format=="str" or format=="text":
3745           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 2980  def mult(arg0,arg1): Line 3798  def mult(arg0,arg1):
3798         """         """
3799         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3800         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3801            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3802         else:         else:
3803            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3804                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3004  class Mult_Symbol(DependendSymbol): Line 3822  class Mult_Symbol(DependendSymbol):
3822         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3823         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3824         """         """
3825         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3826         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3827         if not sh0==sh1:         if not sh0==sh1:
3828            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3829         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 3020  class Mult_Symbol(DependendSymbol): Line 3838  class Mult_Symbol(DependendSymbol):
3838        @type format: C{str}        @type format: C{str}
3839        @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.
3840        @rtype: C{str}        @rtype: C{str}
3841        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3842        """        """
3843        if format=="str" or format=="text":        if format=="str" or format=="text":
3844           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3080  def quotient(arg0,arg1): Line 3898  def quotient(arg0,arg1):
3898         """         """
3899         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3900         if testForZero(args[0]):         if testForZero(args[0]):
3901            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3902         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3903            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3904               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3109  class Quotient_Symbol(DependendSymbol): Line 3927  class Quotient_Symbol(DependendSymbol):
3927         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3928         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3929         """         """
3930         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3931         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3932         if not sh0==sh1:         if not sh0==sh1:
3933            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3934         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 3125  class Quotient_Symbol(DependendSymbol): Line 3943  class Quotient_Symbol(DependendSymbol):
3943        @type format: C{str}        @type format: C{str}
3944        @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.
3945        @rtype: C{str}        @rtype: C{str}
3946        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3947        """        """
3948        if format=="str" or format=="text":        if format=="str" or format=="text":
3949           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3186  def power(arg0,arg1): Line 4004  def power(arg0,arg1):
4004         """         """
4005         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4006         if testForZero(args[0]):         if testForZero(args[0]):
4007            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4008         elif testForZero(args[1]):         elif testForZero(args[1]):
4009            return numarray.ones(args[0],numarray.Float)            return numarray.ones(getShape(args[1]),numarray.Float64)
4010         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4011            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4012         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 3211  class Power_Symbol(DependendSymbol): Line 4029  class Power_Symbol(DependendSymbol):
4029         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4030         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4031         """         """
4032         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4033         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4034         if not sh0==sh1:         if not sh0==sh1:
4035            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4036         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3229  class Power_Symbol(DependendSymbol): Line 4047  class Power_Symbol(DependendSymbol):
4047        @type format: C{str}        @type format: C{str}
4048        @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.
4049        @rtype: C{str}        @rtype: C{str}
4050        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4051        """        """
4052        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4053           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3311  def minimum(*args): Line 4129  def minimum(*args):
4129            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4130      return out      return out
4131    
4132  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4133      """      """
4134      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4135    
4136      @param arg: argument      @param arg: argument
4137      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4138      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4139      @type arg: C{float}      @type minval: C{float} or C{None}
4140      @param maxval: uper range      @param maxval: upper range. If None no upper range is applied
4141      @type arg: C{float}      @type maxval: C{float} or C{None}
4142      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4143               less then maxval are unchanged.               less then maxval are unchanged.
4144      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4145        @raise ValueError: if minval>maxval
4146      """      """
4147      if minval>maxval:      if not minval==None and not maxval==None:
4148         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4149      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4150        if minval == None:
4151            tmp=arg
4152        else:
4153            tmp=maximum(minval,arg)
4154        if maxval == None:
4155            return tmp
4156        else:
4157            return minimum(tmp,maxval)
4158    
4159        
4160  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3344  def inner(arg0,arg1): Line 4171  def inner(arg0,arg1):
4171      @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}
4172      @param arg1: second argument      @param arg1: second argument
4173      @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}
4174      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4175      @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
4176      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4177      """      """
4178      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4179      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4180      if not sh0==sh1:      if not sh0==sh1:
4181          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4182      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4183    
4184    def outer(arg0,arg1):
4185        """
4186        the outer product of the two argument:
4187    
4188        out[t,s]=arg0[t]*arg1[s]
4189    
4190        where
4191    
4192            - s runs through arg0.Shape
4193            - t runs through arg1.Shape
4194    
4195        @param arg0: first argument
4196        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4197        @param arg1: second argument
4198        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4199        @return: the outer product of arg0 and arg1 at each data point
4200        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4201        """
4202        return generalTensorProduct(arg0,arg1,axis_offset=0)
4203    
4204  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4205      """      """
4206        see L{matrix_mult}
4207        """
4208        return matrix_mult(arg0,arg1)
4209    
4210    def matrix_mult(arg0,arg1):
4211        """
4212      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4213    
4214      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4215    
4216            or      or
4217    
4218      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4219    
4220      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.
4221    
4222      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4223      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3374  def matrixmult(arg0,arg1): Line 4227  def matrixmult(arg0,arg1):
4227      @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
4228      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4229      """      """
4230      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4231      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4232      if not len(sh0)==2 :      if not len(sh0)==2 :
4233          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4234      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4235          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4236      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4237    
4238  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4239      """      """
4240      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  
4241      """      """
4242      return generalTensorProduct(arg0,arg1,offset=0)      return tensor_mult(arg0,arg1)
   
4243    
4244  def tensormult(arg0,arg1):  def tensor_mult(arg0,arg1):
4245      """      """
4246      the tensor product of the two argument:      the tensor product of the two argument:
   
4247            
4248      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4249    
4250      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4251    
4252                   or      or
4253    
4254      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4255    
# Line 3419  def tensormult(arg0,arg1): Line 4258  def tensormult(arg0,arg1):
4258    
4259      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]
4260                                
4261                   or      or
4262    
4263      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]
4264    
4265                   or      or
4266    
4267      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]
4268    
4269      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  
4270      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.
4271    
4272      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4273      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3437  def tensormult(arg0,arg1): Line 4276  def tensormult(arg0,arg1):
4276      @return: the tensor product of arg0 and arg1 at each data point      @return: the 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
4278      """      """
4279      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4280      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4281      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4282         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4283      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):
4284         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4285      else:      else:
4286          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4287    
4288  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4289      """      """
4290      generalized tensor product      generalized tensor product
4291    
4292      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4293    
4294      where s runs through arg0.Shape[:arg0.Rank-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[offset:]  
4295    
4296      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]
4297      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4298            - t runs through arg1.Shape[axis_offset:]
4299    
4300      @param arg0: first argument      @param arg0: first argument
4301      @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}
4302      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4303      @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}
4304      @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.
4305      @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 3471  def generalTensorProduct(arg0,arg1,offse Line 4309  def generalTensorProduct(arg0,arg1,offse
4309      # 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
4310      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4311         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4312             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4313         else:         else:
4314             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4315                 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)
4316             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4317             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4318             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4319             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4320             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4321             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4322             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4323             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4324             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4325             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4326             for i0 in range(d0):             for i0 in range(d0):
4327                      for i1 in range(d1):                      for i1 in range(d1):
4328                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4329             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4330             return out             return out
4331      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4332         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4333             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4334         else:         else:
4335             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)
4336      else:            else:      
4337         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4338                                    
4339  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4340     """     """
4341     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4342     """     """
4343     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4344         """         """
4345         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4346    
4347         @param arg0: numerator         @param arg0: first argument
4348         @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}.
4349         @param arg1: denominator         @param arg1: second argument
4350         @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}.
4351         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4352         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4353         """         """
4354         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4355         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4356         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4357         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4358         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4359         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4360         if not sh01==sh10:         if not sh01==sh10:
4361             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)
4362         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])
4363    
4364     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4365        """        """
# Line 3533  class GeneralTensorProduct_Symbol(Depend Line 4371  class GeneralTensorProduct_Symbol(Depend
4371        @type format: C{str}        @type format: C{str}
4372        @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.
4373        @rtype: C{str}        @rtype: C{str}
4374        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4375        """        """
4376        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4377           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])
4378        else:        else:
4379           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)
4380    
# Line 3561  class GeneralTensorProduct_Symbol(Depend Line 4399  class GeneralTensorProduct_Symbol(Depend
4399           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4400           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4401    
4402  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4403      "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!!!"
4404      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4405      shape0=arg0.getShape()[:arg0.getRank()-offset]  
4406      shape01=arg0.getShape()[arg0.getRank()-offset:]  def transposed_matrix_mult(arg0,arg1):
4407      shape10=arg1.getShape()[:offset]      """
4408      shape1=arg1.getShape()[offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4409      if not shape01==shape10:  
4410          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]
4411    
4412      # create return value:      or
4413      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())  
4414      #      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4415      s0=[[]]  
4416      for k in shape0:      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4417            s=[]  
4418            for j in s0:      The first dimension of arg0 and arg1 must match.
4419                  for i in range(k): s.append(j+[slice(i,i)])  
4420            s0=s      @param arg0: first argument of rank 2
4421      s1=[[]]      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4422      for k in shape1:      @param arg1: second argument of at least rank 1
4423            s=[]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4424            for j in s1:      @return: the product of the transposed of arg0 and arg1 at each data point
4425                  for i in range(k): s.append(j+[slice(i,i)])      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4426            s1=s      @raise ValueError: if the shapes of the arguments are not appropriate
4427      s01=[[]]      """
4428      for k in shape01:      sh0=getShape(arg0)
4429            s=[]      sh1=getShape(arg1)
4430            for j in s01:      if not len(sh0)==2 :
4431                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"first argument must have rank 2"
4432            s01=s      if not len(sh1)==2 and not len(sh1)==1:
4433            raise ValueError,"second argument must have rank 1 or 2"
4434      for i0 in s0:      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4435         for i1 in s1:  
4436           s=escript.Scalar(0.,arg0.getFunctionSpace())  def transposed_tensor_mult(arg0,arg1):
4437           for i01 in s01:      """
4438              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      the tensor product of the transposed of the first and the second argument
4439           out.__setitem__(tuple(i0+i1),s)      
4440      return out      for arg0 of rank 2 this is
4441    
4442        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4443    
4444        or
4445    
4446        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4447    
4448      
4449        and for arg0 of rank 4 this is
4450    
4451        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4452                  
4453        or
4454    
4455        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4456    
4457        or
4458    
4459        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4460    
4461        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4462        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4463    
4464        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4465    
4466        @param arg0: first argument of rank 2 or 4
4467        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4468        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4469        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4470        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4471        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4472        """
4473        sh0=getShape(arg0)
4474        sh1=getShape(arg1)
4475        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4476           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4477        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4478           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4479        else:
4480            raise ValueError,"first argument must have rank 2 or 4"
4481    
4482    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4483        """
4484        generalized tensor product of transposed of arg0 and arg1:
4485    
4486        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4487    
4488        where
4489    
4490            - s runs through arg0.Shape[axis_offset:]
4491            - r runs trough arg0.Shape[:axis_offset]
4492            - t runs through arg1.Shape[axis_offset:]
4493    
4494        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4495        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4496    
4497        @param arg0: first argument
4498        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4499        @param arg1: second argument
4500        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4501        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4502        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4503        """
4504        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4505        arg0,arg1=matchType(arg0,arg1)
4506        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4507        if isinstance(arg0,numarray.NumArray):
4508           if isinstance(arg1,Symbol):
4509               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4510           else:
4511               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4512                   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)
4513               arg0_c=arg0.copy()
4514               arg1_c=arg1.copy()
4515               sh0,sh1=arg0.shape,arg1.shape
4516               d0,d1,d01=1,1,1
4517               for i in sh0[axis_offset:]: d0*=i
4518               for i in sh1[axis_offset:]: d1*=i
4519               for i in sh0[:axis_offset]: d01*=i
4520               arg0_c.resize((d01,d0))
4521               arg1_c.resize((d01,d1))
4522               out=numarray.zeros((d0,d1),numarray.Float64)
4523               for i0 in range(d0):
4524                        for i1 in range(d1):
4525                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4526               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4527               return out
4528        elif isinstance(arg0,escript.Data):
4529           if isinstance(arg1,Symbol):
4530               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4531           else:
4532               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4533        else:      
4534           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4535                    
4536    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4537       """
4538       Symbol representing the general tensor product of the transposed of arg0 and arg1
4539       """
4540       def __init__(self,arg0,arg1,axis_offset=0):
4541           """
4542           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4543    
4544           @param arg0: first argument
4545           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4546           @param arg1: second argument
4547           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4548           @raise ValueError: inconsistent dimensions of arguments.
4549           @note: if both arguments have a spatial dimension, they must equal.
4550           """
4551           sh_arg0=getShape(arg0)
4552           sh_arg1=getShape(arg1)
4553           sh01=sh_arg0[:axis_offset]
4554           sh10=sh_arg1[:axis_offset]
4555           sh0=sh_arg0[axis_offset:]
4556           sh1=sh_arg1[axis_offset:]
4557           if not sh01==sh10:
4558               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)
4559           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4560    
4561       def getMyCode(self,argstrs,format="escript"):
4562          """
4563          returns a program code that can be used to evaluate the symbol.
4564    
4565          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4566          @type argstrs: C{list} of length 2 of C{str}.
4567          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4568          @type format: C{str}
4569          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4570          @rtype: C{str}
4571          @raise NotImplementedError: if the requested format is not available
4572          """
4573          if format=="escript" or format=="str" or format=="text":
4574             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4575          else:
4576             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4577    
4578       def substitute(self,argvals):
4579          """
4580          assigns new values to symbols in the definition of the symbol.
4581          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4582    
4583          @param argvals: new values assigned to symbols
4584          @type argvals: C{dict} with keywords of type L{Symbol}.
4585          @return: result of the substitution process. Operations are executed as much as possible.
4586          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4587          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4588          """
4589          if argvals.has_key(self):
4590             arg=argvals[self]
4591             if self.isAppropriateValue(arg):
4592                return arg
4593             else:
4594                raise TypeError,"%s: new value is not appropriate."%str(self)
4595          else:
4596             args=self.getSubstitutedArguments(argvals)
4597             return generalTransposedTensorProduct(args[0],args[1],args[2])
4598    
4599    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4600        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4601        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4602    
4603    def matrix_transposed_mult(arg0,arg1):
4604        """
4605        matrix-transposed(matrix) product of the two argument:
4606    
4607        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4608    
4609        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4610    
4611        The last dimensions of arg0 and arg1 must match.
4612    
4613        @param arg0: first argument of rank 2
4614        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4615        @param arg1: second argument of rank 2
4616        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4617        @return: the product of arg0 and the transposed of arg1 at each data point
4618        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4619        @raise ValueError: if the shapes of the arguments are not appropriate
4620        """
4621        sh0=getShape(arg0)
4622        sh1=getShape(arg1)
4623        if not len(sh0)==2 :
4624            raise ValueError,"first argument must have rank 2"
4625        if not len(sh1)==2 and not len(sh1)==1:
4626            raise ValueError,"second argument must have rank 1 or 2"
4627        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4628    
4629    def tensor_transposed_mult(arg0,arg1):
4630        """
4631        the tensor product of the first and the transpose of the second argument
4632        
4633        for arg0 of rank 2 this is
4634    
4635        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4636    
4637        and for arg0 of rank 4 this is
4638    
4639        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4640                  
4641        or
4642    
4643        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4644    
4645        In the first case the the second dimension of arg0 and arg1 must match and  
4646        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4647    
4648        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4649    
4650        @param arg0: first argument of rank 2 or 4
4651        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4652        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4653        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4654        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4655        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4656        """
4657        sh0=getShape(arg0)
4658        sh1=getShape(arg1)
4659        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4660           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4661        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4662           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4663        else:
4664            raise ValueError,"first argument must have rank 2 or 4"
4665    
4666    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4667        """
4668        generalized tensor product of transposed of arg0 and arg1:
4669    
4670        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4671    
4672        where
4673    
4674            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4675            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4676            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4677    
4678        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4679        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4680    
4681        @param arg0: first argument
4682        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4683        @param arg1: second argument
4684        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4685        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4686        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4687        """
4688        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4689        arg0,arg1=matchType(arg0,arg1)
4690        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4691        if isinstance(arg0,numarray.NumArray):
4692           if isinstance(arg1,Symbol):
4693               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4694           else:
4695               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4696                   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)
4697               arg0_c=arg0.copy()
4698               arg1_c=arg1.copy()
4699               sh0,sh1=arg0.shape,arg1.shape
4700               d0,d1,d01=1,1,1
4701               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4702               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4703               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4704               arg0_c.resize((d0,d01))
4705               arg1_c.resize((d1,d01))
4706               out=numarray.zeros((d0,d1),numarray.Float64)
4707               for i0 in range(d0):
4708                        for i1 in range(d1):
4709                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4710               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4711               return out
4712        elif isinstance(arg0,escript.Data):
4713           if isinstance(arg1,Symbol):
4714               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4715           else:
4716               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4717        else:      
4718           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4719                    
4720    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4721       """
4722       Symbol representing the general tensor product of arg0 and the transpose of arg1
4723       """
4724       def __init__(self,arg0,arg1,axis_offset=0):
4725           """
4726           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4727    
4728           @param arg0: first argument
4729           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4730           @param arg1: second argument
4731           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4732           @raise ValueError: inconsistent dimensions of arguments.
4733           @note: if both arguments have a spatial dimension, they must equal.
4734           """
4735           sh_arg0=getShape(arg0)
4736           sh_arg1=getShape(arg1)
4737           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4738           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4739           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4740           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4741           if not sh01==sh10:
4742               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)
4743           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4744    
4745       def getMyCode(self,argstrs,format="escript"):
4746          """
4747          returns a program code that can be used to evaluate the symbol.
4748    
4749          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4750          @type argstrs: C{list} of length 2 of C{str}.
4751          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4752          @type format: C{str}
4753          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4754          @rtype: C{str}
4755          @raise NotImplementedError: if the requested format is not available
4756          """
4757          if format=="escript" or format=="str" or format=="text":
4758             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4759          else:
4760             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4761    
4762       def substitute(self,argvals):
4763          """
4764          assigns new values to symbols in the definition of the symbol.
4765          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4766    
4767          @param argvals: new values assigned to symbols
4768          @type argvals: C{dict} with keywords of type L{Symbol}.
4769          @return: result of the substitution process. Operations are executed as much as possible.
4770          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4771          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4772          """
4773          if argvals.has_key(self):
4774             arg=argvals[self]
4775             if self.isAppropriateValue(arg):
4776                return arg
4777             else:
4778                raise TypeError,"%s: new value is not appropriate."%str(self)
4779          else:
4780             args=self.getSubstitutedArguments(argvals)
4781             return generalTensorTransposedProduct(args[0],args[1],args[2])
4782    
4783    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4784        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4785        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4786    
4787  #=========================================================  #=========================================================
4788  #   some little helpers  #  functions dealing with spatial dependency
4789  #=========================================================  #=========================================================
4790  def grad(arg,where=None):  def grad(arg,where=None):
4791      """      """
4792      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4793    
4794        If C{g} is the returned object, then
4795    
4796          - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4797          - 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.
4798          - 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.
4799          - 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.
4800    
4801      @param arg:   Data object representing the function which gradient      @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4802                    to be calculated.      @type arg: L{escript.Data} or L{Symbol}
4803      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4804                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4805        @type where: C{None} or L{escript.FunctionSpace}
4806        @return: gradient of arg.
4807        @rtype: L{escript.Data} or L{Symbol}
4808      """      """
4809      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4810         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3621  def grad(arg,where=None): Line 4814  def grad(arg,where=None):
4814         else:         else:
4815            return arg._grad(where)            return arg._grad(where)
4816      else:      else:
4817        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4818    
4819    class Grad_Symbol(DependendSymbol):
4820       """
4821       L{Symbol} representing the result of the gradient operator
4822       """
4823       def __init__(self,arg,where=None):
4824          """
4825          initialization of gradient L{Symbol} with argument arg
4826          @param arg: argument of function
4827          @type arg: L{Symbol}.
4828          @param where: FunctionSpace in which the gradient will be calculated.
4829                      If not present or C{None} an appropriate default is used.
4830          @type where: C{None} or L{escript.FunctionSpace}
4831          """
4832          d=arg.getDim()
4833          if d==None:
4834             raise ValueError,"argument must have a spatial dimension"
4835          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4836    
4837       def getMyCode(self,argstrs,format="escript"):
4838          """
4839          returns a program code that can be used to evaluate the symbol.
4840    
4841          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4842          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4843          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4844          @type format: C{str}
4845          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4846          @rtype: C{str}
4847          @raise NotImplementedError: if the requested format is not available
4848          """
4849          if format=="escript" or format=="str"  or format=="text":
4850             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4851          else:
4852             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4853    
4854       def substitute(self,argvals):
4855          """
4856          assigns new values to symbols in the definition of the symbol.
4857          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4858    
4859          @param argvals: new values assigned to symbols
4860          @type argvals: C{dict} with keywords of type L{Symbol}.
4861          @return: result of the substitution process. Operations are executed as much as possible.
4862          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4863          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4864          """
4865          if argvals.has_key(self):
4866             arg=argvals[self]
4867             if self.isAppropriateValue(arg):
4868                return arg
4869             else:
4870                raise TypeError,"%s: new value is not appropriate."%str(self)
4871          else:
4872             arg=self.getSubstitutedArguments(argvals)
4873             return grad(arg[0],where=arg[1])
4874    
4875       def diff(self,arg):
4876          """
4877          differential of this object
4878    
4879          @param arg: the derivative is calculated with respect to arg
4880          @type arg: L{escript.Symbol}
4881          @return: derivative with respect to C{arg}
4882          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4883          """
4884          if arg==self:
4885             return identity(self.getShape())
4886          else:
4887             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4888    
4889  def integrate(arg,where=None):  def integrate(arg,where=None):
4890      """      """
4891      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}
4892      its domain.      before integration.
4893    
4894      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4895        @type arg: L{escript.Data} or L{Symbol}
4896      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4897                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4898        @type where: C{None} or L{escript.FunctionSpace}
4899        @return: integral of arg.
4900        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4901      """      """
4902      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):  
4903         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4904      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4905         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
# Line 3658  def integrate(arg,where=None): Line 4910  def integrate(arg,where=None):
4910      else:      else:
4911        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4912    
4913    class Integrate_Symbol(DependendSymbol):
4914       """
4915       L{Symbol} representing the result of the spatial integration operator
4916       """
4917       def __init__(self,arg,where=None):
4918          """
4919          initialization of integration L{Symbol} with argument arg
4920          @param arg: argument of the integration
4921          @type arg: L{Symbol}.
4922          @param where: FunctionSpace in which the integration will be calculated.
4923                      If not present or C{None} an appropriate default is used.
4924          @type where: C{None} or L{escript.FunctionSpace}
4925          """
4926          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4927    
4928       def getMyCode(self,argstrs,format="escript"):
4929          """
4930          returns a program code that can be used to evaluate the symbol.
4931    
4932          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4933          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4934          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4935          @type format: C{str}
4936          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4937          @rtype: C{str}
4938          @raise NotImplementedError: if the requested format is not available
4939          """
4940          if format=="escript" or format=="str"  or format=="text":
4941             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4942          else:
4943             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4944    
4945       def substitute(self,argvals):
4946          """
4947          assigns new values to symbols in the definition of the symbol.
4948          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4949    
4950          @param argvals: new values assigned to symbols
4951          @type argvals: C{dict} with keywords of type L{Symbol}.
4952          @return: result of the substitution process. Operations are executed as much as possible.
4953          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4954          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4955          """
4956          if argvals.has_key(self):
4957             arg=argvals[self]
4958             if self.isAppropriateValue(arg):
4959                return arg
4960             else:
4961                raise TypeError,"%s: new value is not appropriate."%str(self)
4962          else:
4963             arg=self.getSubstitutedArguments(argvals)
4964             return integrate(arg[0],where=arg[1])
4965    
4966       def diff(self,arg):
4967          """
4968          differential of this object
4969    
4970          @param arg: the derivative is calculated with respect to arg
4971          @type arg: L{escript.Symbol}
4972          @return: derivative with respect to C{arg}
4973          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4974          """
4975          if arg==self:
4976             return identity(self.getShape())
4977          else:
4978             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4979    
4980    
4981  def interpolate(arg,where):  def interpolate(arg,where):
4982      """      """
4983      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
4984    
4985      @param arg:    interpolant      @param arg: interpolant
4986      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
4987        @param where: FunctionSpace to be interpolated to
4988        @type where: L{escript.FunctionSpace}
4989        @return: interpolated argument
4990        @rtype: C{escript.Data} or L{Symbol}
4991      """      """
4992      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4993         return Interpolated_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
4994      else:      else:
4995         return escript.Data(arg,where)         return escript.Data(arg,where)
4996    
4997  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
4998      """     """
4999      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
5000       """
5001       def __init__(self,arg,where):
5002          """
5003          initialization of interpolation L{Symbol} with argument arg
5004          @param arg: argument of the interpolation
5005          @type arg: L{Symbol}.
5006          @param where: FunctionSpace into which the argument is interpolated.
5007          @type where: L{escript.FunctionSpace}
5008          """
5009          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
5010    
5011      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
5012                    be calculated.        """
5013      @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)  
5014    
5015  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
5016      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
5017      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.
5018          @type format: C{str}
5019          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
5020          @rtype: C{str}
5021          @raise NotImplementedError: if the requested format is not available
5022          """
5023          if format=="escript" or format=="str"  or format=="text":
5024             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
5025          else:
5026             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
5027    
5028      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
5029                    to be calculated.        """
5030      """        assigns new values to symbols in the definition of the symbol.
5031      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))  
5032    
5033  #=============================        @param argvals: new values assigned to symbols
5034  #        @type argvals: C{dict} with keywords of type L{Symbol}.
5035  # 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.
5036  # 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
5037  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
5038          """
5039          if argvals.has_key(self):
5040             arg=argvals[self]
5041             if self.isAppropriateValue(arg):
5042                return arg
5043             else:
5044                raise TypeError,"%s: new value is not appropriate."%str(self)
5045          else:
5046             arg=self.getSubstitutedArguments(argvals)
5047             return interpolate(arg[0],where=arg[1])
5048    
5049  # functions involving the underlying Domain:     def diff(self,arg):
5050          """
5051          differential of this object
5052    
5053          @param arg: the derivative is calculated with respect to arg
5054          @type arg: L{escript.Symbol}
5055          @return: derivative with respect to C{arg}
5056          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
5057          """
5058          if arg==self:
5059             return identity(self.getShape())
5060          else:
5061             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5062    
5063  def transpose(arg,axis=None):  
5064    def div(arg,where=None):
5065      """      """
5066      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
5067    
5068      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
5069        @type arg: L{escript.Data} or L{Symbol}
5070        @param where: FunctionSpace in which the divergence will be calculated.
5071                      If not present or C{None} an appropriate default is used.
5072        @type where: C{None} or L{escript.FunctionSpace}
5073        @return: divergence of arg.
5074        @rtype: L{escript.Data} or L{Symbol}
5075      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
5076      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5077         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
5078      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
5079         # 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)  
5080      else:      else:
5081         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
5082        if not arg.getShape()==(dim,):
5083          raise ValueError,"div: expected shape is (%s,)"%dim
5084        return trace(grad(arg,where))
5085    
5086  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
5087      """      """
5088      Return      returns the jump of arg across the continuity of the domain
5089    
5090      @param arg:      @param arg: argument
5091        @type arg: L{escript.Data} or L{Symbol}
5092        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5093                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5094        @type domain: C{None} or L{escript.Domain}
5095        @return: jump of arg
5096        @rtype: L{escript.Data} or L{Symbol}
5097      """      """
5098      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
5099         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)  
5100    
5101    def L2(arg):
5102        """
5103        returns the L2 norm of arg at where
5104        
5105        @param arg: function which L2 to be calculated.
5106        @type arg: L{escript.Data} or L{Symbol}
5107        @return: L2 norm of arg.
5108        @rtype: L{float} or L{Symbol}
5109        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5110        """
5111        return sqrt(integrate(inner(arg,arg)))
5112    #=============================
5113    #
5114    
5115  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5116      """      """
5117      resorts the component of arg according to index      resorts the component of arg according to index
5118    
5119      """      """
5120      pass      raise NotImplementedError
5121  #  #
5122  # $Log: util.py,v $  # $Log: util.py,v $
5123  # Revision 1.14.2.16  2005/10/19 06:09:57  gross  # Revision 1.14.2.16  2005/10/19 06:09:57  gross

Legend:
Removed from v.396  
changed lines
  Added in v.1044

  ViewVC Help
Powered by ViewVC 1.1.26