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

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

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

revision 341 by gross, Mon Dec 12 05:26:10 2005 UTC revision 912 by gross, Wed Dec 6 03:29:49 2006 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2    
3  """  """
4  Utility functions for escript  Utility functions for escript
5    
 @remark:  This module is under construction and is still tested!!!  
   
6  @var __author__: name of author  @var __author__: name of author
7  @var __licence__: licence agreement  @var __copyright__: copyrights
8    @var __license__: licence agreement
9  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
10  @var __version__: version  @var __version__: version
11  @var __date__: date of the version  @var __date__: date of the version
12  """  """
13                                                                                                                                                                                                                                                                                                                                                                                                            
14  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
15  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                        http://www.access.edu.au
17                    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19                 http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
21  __version__="$Revision: 329 $"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
23    
24    
# Line 32  import math Line 26  import math
26  import numarray  import numarray
27  import escript  import escript
28  import os  import os
29    from esys.escript import C_GeneralTensorProduct
 # missing tests:  
   
 # def pokeShape(arg):  
 # def pokeDim(arg):  
 # def commonShape(arg0,arg1):  
 # def commonDim(*args):  
 # def testForZero(arg):  
 # def matchType(arg0=0.,arg1=0.):  
 # def matchShape(arg0,arg1):  
   
 # def maximum(arg0,arg1):  
 # def minimum(arg0,arg1):  
   
 # def transpose(arg,axis=None):  
 # def trace(arg,axis0=0,axis1=1):  
 # def reorderComponents(arg,index):  
   
 # def integrate(arg,where=None):  
 # def interpolate(arg,where):  
 # def div(arg,where=None):  
 # def grad(arg,where=None):  
   
 #  
 # slicing: get  
 #          set  
 #  
 # and derivatives  
30    
31  #=========================================================  #=========================================================
32  #   some helpers:  #   some helpers:
# Line 68  def saveVTK(filename,domain=None,**data) Line 35  def saveVTK(filename,domain=None,**data)
35      """      """
36      writes a L{Data} objects into a files using the the VTK XML file format.      writes a L{Data} objects into a files using the the VTK XML file format.
37    
38      Example:      Example::
39    
40         tmp=Scalar(..)         tmp=Scalar(..)
41         v=Vector(..)         v=Vector(..)
# Line 96  def saveDX(filename,domain=None,**data): Line 63  def saveDX(filename,domain=None,**data):
63      """      """
64      writes a L{Data} objects into a files using the the DX file format.      writes a L{Data} objects into a files using the the DX file format.
65    
66      Example:      Example::
67    
68         tmp=Scalar(..)         tmp=Scalar(..)
69         v=Vector(..)         v=Vector(..)
# Line 125  def kronecker(d=3): Line 92  def kronecker(d=3):
92     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
93    
94     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
95     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
96     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
97     @rtype d: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2.
    @remark: the function is identical L{identity}  
98     """     """
99     return identityTensor(d)     return identityTensor(d)
100    
# Line 143  def identity(shape=()): Line 109  def identity(shape=()):
109     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
110     """     """
111     if len(shape)>0:     if len(shape)>0:
112        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
113        if len(shape)==1:        if len(shape)==1:
114            for i0 in range(shape[0]):            for i0 in range(shape[0]):
115               out[i0,i0]=1.               out[i0,i0]=1.
   
116        elif len(shape)==2:        elif len(shape)==2:
117            for i0 in range(shape[0]):            for i0 in range(shape[0]):
118               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 163  def identityTensor(d=3): Line 128  def identityTensor(d=3):
128     return the dxd identity matrix     return the dxd identity matrix
129    
130     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
131     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
132     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
133     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
134     """     """
135     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
136        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
137     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
138           return identity((d.getDim(),))
139       else:
140           return identity((d,))
141    
142  def identityTensor4(d=3):  def identityTensor4(d=3):
143     """     """
# Line 178  def identityTensor4(d=3): Line 146  def identityTensor4(d=3):
146     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
147     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
148     @return: the object u of rank 4 with M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise     @return: the object u of rank 4 with M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise
149     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
150     """     """
151     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
152        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
153     return identity((d,d))     elif isinstance(d,escript.Domain):
154           return identity((d.getDim(),d.getDim()))
155       else:
156           return identity((d,d))
157    
158  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
159     """     """
# Line 191  def unitVector(i=0,d=3): Line 162  def unitVector(i=0,d=3):
162     @param i: index     @param i: index
163     @type i: C{int}     @type i: C{int}
164     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
165     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
166     @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise     @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise
167     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
168     """     """
169     return kronecker(d)[i]     return kronecker(d)[i]
170    
# Line 249  def inf(arg): Line 220  def inf(arg):
220    
221      @param arg: argument      @param arg: argument
222      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
223      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
224      @rtype: C{float}      @rtype: C{float}
225      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
226      """      """
# Line 268  def inf(arg): Line 239  def inf(arg):
239  #=========================================================================  #=========================================================================
240  #   some little helpers  #   some little helpers
241  #=========================================================================  #=========================================================================
242  def pokeShape(arg):  def getRank(arg):
243        """
244        identifies the rank of its argument
245    
246        @param arg: a given object
247        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
248        @return: the rank of the argument
249        @rtype: C{int}
250        @raise TypeError: if type of arg cannot be processed
251        """
252    
253        if isinstance(arg,numarray.NumArray):
254            return arg.rank
255        elif isinstance(arg,escript.Data):
256            return arg.getRank()
257        elif isinstance(arg,float):
258            return 0
259        elif isinstance(arg,int):
260            return 0
261        elif isinstance(arg,Symbol):
262            return arg.getRank()
263        else:
264          raise TypeError,"getShape: cannot identify shape"
265    def getShape(arg):
266      """      """
267      identifies the shape of its argument      identifies the shape of its argument
268    
# Line 290  def pokeShape(arg): Line 284  def pokeShape(arg):
284      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
285          return arg.getShape()          return arg.getShape()
286      else:      else:
287        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
288    
289  def pokeDim(arg):  def pokeDim(arg):
290      """      """
# Line 313  def commonShape(arg0,arg1): Line 307  def commonShape(arg0,arg1):
307      """      """
308      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.
309    
310      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
311      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
312      @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.
313      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
314      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
315      """      """
316      sh0=pokeShape(arg0)      sh0=getShape(arg0)
317      sh1=pokeShape(arg1)      sh1=getShape(arg1)
318      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
319         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
320               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
# Line 338  def commonDim(*args): Line 332  def commonDim(*args):
332      """      """
333      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.
334    
335      @param *args: given objects      @param args: given objects
336      @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
337               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
338      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 360  def testForZero(arg): Line 354  def testForZero(arg):
354    
355      @param arg: a given object      @param arg: a given object
356      @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}
357      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
358      @rtype : C{bool}      @rtype: C{bool}
359      """      """
360      try:      if isinstance(arg,numarray.NumArray):
361           return not Lsup(arg)>0.
362        elif isinstance(arg,escript.Data):
363           return False
364        elif isinstance(arg,float):
365           return not Lsup(arg)>0.
366        elif isinstance(arg,int):
367         return not Lsup(arg)>0.         return not Lsup(arg)>0.
368      except TypeError:      elif isinstance(arg,Symbol):
369           return False
370        else:
371         return False         return False
372    
373  def matchType(arg0=0.,arg1=0.):  def matchType(arg0=0.,arg1=0.):
# Line 386  def matchType(arg0=0.,arg1=0.): Line 388  def matchType(arg0=0.,arg1=0.):
388         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
389            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
390         elif isinstance(arg1,float):         elif isinstance(arg1,float):
391            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
392         elif isinstance(arg1,int):         elif isinstance(arg1,int):
393            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
394         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
395            pass            pass
396         else:         else:
# Line 412  def matchType(arg0=0.,arg1=0.): Line 414  def matchType(arg0=0.,arg1=0.):
414         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
415            pass            pass
416         elif isinstance(arg1,float):         elif isinstance(arg1,float):
417            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
418         elif isinstance(arg1,int):         elif isinstance(arg1,int):
419            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
420         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
421            pass            pass
422         else:         else:
423            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
424      elif isinstance(arg0,float):      elif isinstance(arg0,float):
425         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
426            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
427         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
428            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
429         elif isinstance(arg1,float):         elif isinstance(arg1,float):
430            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
431            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
432         elif isinstance(arg1,int):         elif isinstance(arg1,int):
433            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
434            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
435         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
436            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
437         else:         else:
438            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
439      elif isinstance(arg0,int):      elif isinstance(arg0,int):
440         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
441            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
442         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
443            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
444         elif isinstance(arg1,float):         elif isinstance(arg1,float):
445            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
446            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
447         elif isinstance(arg1,int):         elif isinstance(arg1,int):
448            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
449            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
450         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
451            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
452         else:         else:
453            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
454      else:      else:
# Line 456  def matchType(arg0=0.,arg1=0.): Line 458  def matchType(arg0=0.,arg1=0.):
458    
459  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
460      """      """
461            return representations of arg0 amd arg1 which ahve the same shape
462    
463      If shape is not given the shape "largest" shape of args is used.      @param arg0: a given object
464        @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
465      @param args: a given ob      @param arg1: a given object
466      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
467      @return: True if the argument is identical to zero.      @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
468      @rtype: C{list} of C{int}      @rtype: C{tuple}
469      """      """
470      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
471      sh0=pokeShape(arg0)      sh0=getShape(arg0)
472      sh1=pokeShape(arg1)      sh1=getShape(arg1)
473      if len(sh0)<len(sh):      if len(sh0)<len(sh):
474         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
475      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
476         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
477      else:      else:
478         return arg0,arg1         return arg0,arg1
479  #=========================================================  #=========================================================
# Line 491  class Symbol(object): Line 493  class Symbol(object):
493         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
494         symbols or any other object.         symbols or any other object.
495    
496         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
497         @type arg: C{list}         @type args: C{list}
498         @param shape: the shape         @param shape: the shape
499         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
500         @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
# Line 535  class Symbol(object): Line 537  class Symbol(object):
537         """         """
538         the shape of the symbol.         the shape of the symbol.
539    
540         @return : the shape of the symbol.         @return: the shape of the symbol.
541         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
542         """         """
543         return self.__shape         return self.__shape
# Line 544  class Symbol(object): Line 546  class Symbol(object):
546         """         """
547         the spatial dimension         the spatial dimension
548    
549         @return : the spatial dimension         @return: the spatial dimension
550         @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.
551         """         """
552         return self.__dim         return self.__dim
# Line 568  class Symbol(object): Line 570  class Symbol(object):
570         """         """
571         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.
572    
573         @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].
574         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
575         @rtype: C{list} of objects         @rtype: C{list} of objects
576         @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.
577         """         """
578         out=[]         out=[]
579         for a in self.getArgument():         for a in self.getArgument():
# Line 595  class Symbol(object): Line 597  class Symbol(object):
597            if isinstance(a,Symbol):            if isinstance(a,Symbol):
598               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
599            else:            else:
600                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
601                if len(s)>0:                if len(s)>0:
602                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
603                else:                else:
604                   out.append(a)                   out.append(a)
605         return out         return out
# Line 687  class Symbol(object): Line 689  class Symbol(object):
689         else:         else:
690            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
691            if len(s)>0:            if len(s)>0:
692               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
693            else:            else:
694               return 0.               return 0.
695    
# Line 695  class Symbol(object): Line 697  class Symbol(object):
697         """         """
698         returns -self.         returns -self.
699    
700         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
701         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
702         """         """
703         return self*(-1.)         return self*(-1.)
# Line 704  class Symbol(object): Line 706  class Symbol(object):
706         """         """
707         returns +self.         returns +self.
708    
709         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
710         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
711         """         """
712         return self*(1.)         return self*(1.)
713    
714     def __abs__(self):     def __abs__(self):
715         """         """
716         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
717         """         """
718         return Abs_Symbol(self)         return Abs_Symbol(self)
719    
# Line 721  class Symbol(object): Line 723  class Symbol(object):
723    
724         @param other: object to be added to this object         @param other: object to be added to this object
725         @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}.
726         @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}
727         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
728         """         """
729         return add(self,other)         return add(self,other)
# Line 732  class Symbol(object): Line 734  class Symbol(object):
734    
735         @param other: object this object is added to         @param other: object this object is added to
736         @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}.
737         @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
738         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
739         """         """
740         return add(other,self)         return add(other,self)
# Line 743  class Symbol(object): Line 745  class Symbol(object):
745    
746         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
747         @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}.
748         @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
749         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
750         """         """
751         return add(self,-other)         return add(self,-other)
# Line 754  class Symbol(object): Line 756  class Symbol(object):
756    
757         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
758         @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}.
759         @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}.
760         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
761         """         """
762         return add(-self,other)         return add(-self,other)
# Line 765  class Symbol(object): Line 767  class Symbol(object):
767    
768         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
769         @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}.
770         @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}.
771         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
772         """         """
773         return mult(self,other)         return mult(self,other)
# Line 776  class Symbol(object): Line 778  class Symbol(object):
778    
779         @param other: object this object is multiplied with         @param other: object this object is multiplied with
780         @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}.
781         @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.
782         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
783         """         """
784         return mult(other,self)         return mult(other,self)
# Line 787  class Symbol(object): Line 789  class Symbol(object):
789    
790         @param other: object dividing this object         @param other: object dividing this object
791         @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}.
792         @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}
793         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
794         """         """
795         return quotient(self,other)         return quotient(self,other)
# Line 798  class Symbol(object): Line 800  class Symbol(object):
800    
801         @param other: object dividing this object         @param other: object dividing this object
802         @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}.
803         @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
804         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
805         """         """
806         return quotient(other,self)         return quotient(other,self)
# Line 809  class Symbol(object): Line 811  class Symbol(object):
811    
812         @param other: exponent         @param other: exponent
813         @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}.
814         @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}
815         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
816         """         """
817         return power(self,other)         return power(self,other)
# Line 820  class Symbol(object): Line 822  class Symbol(object):
822    
823         @param other: basis         @param other: basis
824         @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}.
825         @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
826         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
827         """         """
828         return power(other,self)         return power(other,self)
829    
830       def __getitem__(self,index):
831           """
832           returns the slice defined by index
833    
834           @param index: defines a
835           @type index: C{slice} or C{int} or a C{tuple} of them
836           @return: a L{Symbol} representing the slice defined by index
837           @rtype: L{DependendSymbol}
838           """
839           return GetSlice_Symbol(self,index)
840    
841  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
842     """     """
843     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.
844     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  
845        
846     Example:     Example::
847        
848     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
849     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
850     print u1==u2       print u1==u2
851     False       False
852        
853        but       but::
854    
855     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
856     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
857     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
858     print u1==u2, u1==u3       print u1==u2, u1==u3
859     True False       True False
860    
861     @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.
862     """     """
# Line 875  class DependendSymbol(Symbol): Line 888  class DependendSymbol(Symbol):
888  #=========================================================  #=========================================================
889  #  Unary operations prserving the shape  #  Unary operations prserving the shape
890  #========================================================  #========================================================
891    class GetSlice_Symbol(DependendSymbol):
892       """
893       L{Symbol} representing getting a slice for a L{Symbol}
894       """
895       def __init__(self,arg,index):
896          """
897          initialization of wherePositive L{Symbol} with argument arg
898          @param arg: argument
899          @type arg: L{Symbol}.
900          @param index: defines index
901          @type index: C{slice} or C{int} or a C{tuple} of them
902          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
903          @raises ValueError: if a step is given
904          """
905          if not isinstance(index,tuple): index=(index,)
906          if len(index)>arg.getRank():
907               raise IndexError,"GetSlice_Symbol: index out of range."
908          sh=()
909          index2=()
910          for i in range(len(index)):
911             ix=index[i]
912             if isinstance(ix,int):
913                if ix<0 or ix>=arg.getShape()[i]:
914                   raise ValueError,"GetSlice_Symbol: index out of range."
915                index2=index2+(ix,)
916             else:
917               if not ix.step==None:
918                 raise ValueError,"GetSlice_Symbol: steping is not supported."
919               if ix.start==None:
920                  s=0
921               else:
922                  s=ix.start
923               if ix.stop==None:
924                  e=arg.getShape()[i]
925               else:
926                  e=ix.stop
927                  if e>arg.getShape()[i]:
928                     raise IndexError,"GetSlice_Symbol: index out of range."
929               index2=index2+(slice(s,e),)
930               if e>s:
931                   sh=sh+(e-s,)
932               elif s>e:
933                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
934          for i in range(len(index),arg.getRank()):
935              index2=index2+(slice(0,arg.getShape()[i]),)
936              sh=sh+(arg.getShape()[i],)
937          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
938    
939       def getMyCode(self,argstrs,format="escript"):
940          """
941          returns a program code that can be used to evaluate the symbol.
942    
943          @param argstrs: gives for each argument a string representing the argument for the evaluation.
944          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
945          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
946          @type format: C{str}
947          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
948          @rtype: C{str}
949          @raise NotImplementedError: if the requested format is not available
950          """
951          if format=="escript" or format=="str"  or format=="text":
952             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
953          else:
954             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
955    
956       def substitute(self,argvals):
957          """
958          assigns new values to symbols in the definition of the symbol.
959          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
960    
961          @param argvals: new values assigned to symbols
962          @type argvals: C{dict} with keywords of type L{Symbol}.
963          @return: result of the substitution process. Operations are executed as much as possible.
964          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
965          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
966          """
967          if argvals.has_key(self):
968             arg=argvals[self]
969             if self.isAppropriateValue(arg):
970                return arg
971             else:
972                raise TypeError,"%s: new value is not appropriate."%str(self)
973          else:
974             args=self.getSubstitutedArguments(argvals)
975             arg=args[0]
976             index=args[1]
977             return arg.__getitem__(index)
978    
979  def log10(arg):  def log10(arg):
980     """     """
981     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
982    
983     @param arg: argument     @param arg: argument
984     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
985     @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.
986     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
987     """     """
988     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 903  def wherePositive(arg): Line 1004  def wherePositive(arg):
1004    
1005     @param arg: argument     @param arg: argument
1006     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1007     @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.
1008     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1009     """     """
1010     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1011        if arg.rank==0:        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1012           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1013             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))  
1014     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1015        return arg._wherePositive()        return arg._wherePositive()
1016     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 953  class WherePositive_Symbol(DependendSymb Line 1050  class WherePositive_Symbol(DependendSymb
1050        @type format: C{str}        @type format: C{str}
1051        @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.
1052        @rtype: C{str}        @rtype: C{str}
1053        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1054        """        """
1055        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1056            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 989  def whereNegative(arg): Line 1086  def whereNegative(arg):
1086    
1087     @param arg: argument     @param arg: argument
1088     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1089     @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.
1090     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1091     """     """
1092     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1093        if arg.rank==0:        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1094           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1095             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))  
1096     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1097        return arg._whereNegative()        return arg._whereNegative()
1098     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1039  class WhereNegative_Symbol(DependendSymb Line 1132  class WhereNegative_Symbol(DependendSymb
1132        @type format: C{str}        @type format: C{str}
1133        @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.
1134        @rtype: C{str}        @rtype: C{str}
1135        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1136        """        """
1137        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1138            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1075  def whereNonNegative(arg): Line 1168  def whereNonNegative(arg):
1168    
1169     @param arg: argument     @param arg: argument
1170     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1171     @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.
1172     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1173     """     """
1174     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1175        if arg.rank==0:        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1176           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1177             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))  
1178     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1179        return arg._whereNonNegative()        return arg._whereNonNegative()
1180     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1109  def whereNonPositive(arg): Line 1198  def whereNonPositive(arg):
1198    
1199     @param arg: argument     @param arg: argument
1200     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1201     @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.
1202     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1203     """     """
1204     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1205        if arg.rank==0:        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1206           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1207             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.  
1208     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1209        return arg._whereNonPositive()        return arg._whereNonPositive()
1210     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1145  def whereZero(arg,tol=0.): Line 1230  def whereZero(arg,tol=0.):
1230     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1231     @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.
1232     @type tol: C{float}     @type tol: C{float}
1233     @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.
1234     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1235     """     """
1236     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1237        if arg.rank==0:        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1238           if abs(arg)<=tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1239             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
1240     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1241        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1242     elif isinstance(arg,float):     elif isinstance(arg,float):
1243        if abs(arg)<=tol:        if abs(arg)<=tol:
1244          return 1.          return 1.
# Line 1198  class WhereZero_Symbol(DependendSymbol): Line 1276  class WhereZero_Symbol(DependendSymbol):
1276        @type format: C{str}        @type format: C{str}
1277        @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.
1278        @rtype: C{str}        @rtype: C{str}
1279        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1280        """        """
1281        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1282           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1232  def whereNonZero(arg,tol=0.): Line 1310  def whereNonZero(arg,tol=0.):
1310    
1311     @param arg: argument     @param arg: argument
1312     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1313     @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.
1314     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1315     """     """
1316     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1317        if arg.rank==0:        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1318          if abs(arg)>tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1319             return numarray.array(1.)        return out
         else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
1320     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1321        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1322     elif isinstance(arg,float):     elif isinstance(arg,float):
1323        if abs(arg)>tol:        if abs(arg)>tol:
1324          return 1.          return 1.
# Line 1263  def whereNonZero(arg,tol=0.): Line 1334  def whereNonZero(arg,tol=0.):
1334     else:     else:
1335        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1336    
1337    def erf(arg):
1338       """
1339       returns erf of argument arg
1340    
1341       @param arg: argument
1342       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1343       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1344       @raises TypeError: if the type of the argument is not expected.
1345       """
1346       if isinstance(arg,escript.Data):
1347          return arg._erf()
1348       else:
1349          raise TypeError,"erf: Unknown argument type."
1350    
1351  def sin(arg):  def sin(arg):
1352     """     """
1353     returns sine of argument arg     returns sine of argument arg
1354    
1355     @param arg: argument     @param arg: argument
1356     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1357     @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.
1358     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1359     """     """
1360     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1307  class Sin_Symbol(DependendSymbol): Line 1392  class Sin_Symbol(DependendSymbol):
1392        @type format: C{str}        @type format: C{str}
1393        @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.
1394        @rtype: C{str}        @rtype: C{str}
1395        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1396        """        """
1397        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1398            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1359  def cos(arg): Line 1444  def cos(arg):
1444    
1445     @param arg: argument     @param arg: argument
1446     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1447     @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.
1448     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1449     """     """
1450     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1397  class Cos_Symbol(DependendSymbol): Line 1482  class Cos_Symbol(DependendSymbol):
1482        @type format: C{str}        @type format: C{str}
1483        @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.
1484        @rtype: C{str}        @rtype: C{str}
1485        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1486        """        """
1487        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1488            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1449  def tan(arg): Line 1534  def tan(arg):
1534    
1535     @param arg: argument     @param arg: argument
1536     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1537     @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.
1538     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1539     """     """
1540     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1487  class Tan_Symbol(DependendSymbol): Line 1572  class Tan_Symbol(DependendSymbol):
1572        @type format: C{str}        @type format: C{str}
1573        @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.
1574        @rtype: C{str}        @rtype: C{str}
1575        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1576        """        """
1577        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1578            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1539  def asin(arg): Line 1624  def asin(arg):
1624    
1625     @param arg: argument     @param arg: argument
1626     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1627     @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.
1628     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1629     """     """
1630     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1577  class Asin_Symbol(DependendSymbol): Line 1662  class Asin_Symbol(DependendSymbol):
1662        @type format: C{str}        @type format: C{str}
1663        @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.
1664        @rtype: C{str}        @rtype: C{str}
1665        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1666        """        """
1667        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1668            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1629  def acos(arg): Line 1714  def acos(arg):
1714    
1715     @param arg: argument     @param arg: argument
1716     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1717     @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.
1718     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1719     """     """
1720     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1667  class Acos_Symbol(DependendSymbol): Line 1752  class Acos_Symbol(DependendSymbol):
1752        @type format: C{str}        @type format: C{str}
1753        @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.
1754        @rtype: C{str}        @rtype: C{str}
1755        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1756        """        """
1757        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1758            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1719  def atan(arg): Line 1804  def atan(arg):
1804    
1805     @param arg: argument     @param arg: argument
1806     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1807     @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.
1808     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1809     """     """
1810     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1757  class Atan_Symbol(DependendSymbol): Line 1842  class Atan_Symbol(DependendSymbol):
1842        @type format: C{str}        @type format: C{str}
1843        @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.
1844        @rtype: C{str}        @rtype: C{str}
1845        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1846        """        """
1847        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1848            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1809  def sinh(arg): Line 1894  def sinh(arg):
1894    
1895     @param arg: argument     @param arg: argument
1896     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1897     @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.
1898     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1899     """     """
1900     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1847  class Sinh_Symbol(DependendSymbol): Line 1932  class Sinh_Symbol(DependendSymbol):
1932        @type format: C{str}        @type format: C{str}
1933        @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.
1934        @rtype: C{str}        @rtype: C{str}
1935        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1936        """        """
1937        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1938            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1899  def cosh(arg): Line 1984  def cosh(arg):
1984    
1985     @param arg: argument     @param arg: argument
1986     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1987     @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.
1988     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1989     """     """
1990     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1937  class Cosh_Symbol(DependendSymbol): Line 2022  class Cosh_Symbol(DependendSymbol):
2022        @type format: C{str}        @type format: C{str}
2023        @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.
2024        @rtype: C{str}        @rtype: C{str}
2025        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2026        """        """
2027        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2028            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1989  def tanh(arg): Line 2074  def tanh(arg):
2074    
2075     @param arg: argument     @param arg: argument
2076     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2077     @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.
2078     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2079     """     """
2080     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2027  class Tanh_Symbol(DependendSymbol): Line 2112  class Tanh_Symbol(DependendSymbol):
2112        @type format: C{str}        @type format: C{str}
2113        @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.
2114        @rtype: C{str}        @rtype: C{str}
2115        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2116        """        """
2117        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2118            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2079  def asinh(arg): Line 2164  def asinh(arg):
2164    
2165     @param arg: argument     @param arg: argument
2166     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2167     @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.
2168     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2169     """     """
2170     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2117  class Asinh_Symbol(DependendSymbol): Line 2202  class Asinh_Symbol(DependendSymbol):
2202        @type format: C{str}        @type format: C{str}
2203        @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.
2204        @rtype: C{str}        @rtype: C{str}
2205        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2206        """        """
2207        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2208            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2169  def acosh(arg): Line 2254  def acosh(arg):
2254    
2255     @param arg: argument     @param arg: argument
2256     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2257     @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.
2258     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2259     """     """
2260     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2207  class Acosh_Symbol(DependendSymbol): Line 2292  class Acosh_Symbol(DependendSymbol):
2292        @type format: C{str}        @type format: C{str}
2293        @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.
2294        @rtype: C{str}        @rtype: C{str}
2295        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2296        """        """
2297        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2298            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2259  def atanh(arg): Line 2344  def atanh(arg):
2344    
2345     @param arg: argument     @param arg: argument
2346     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2347     @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.
2348     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2349     """     """
2350     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2297  class Atanh_Symbol(DependendSymbol): Line 2382  class Atanh_Symbol(DependendSymbol):
2382        @type format: C{str}        @type format: C{str}
2383        @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.
2384        @rtype: C{str}        @rtype: C{str}
2385        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2386        """        """
2387        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2388            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2349  def exp(arg): Line 2434  def exp(arg):
2434    
2435     @param arg: argument     @param arg: argument
2436     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2437     @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.
2438     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2439     """     """
2440     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2387  class Exp_Symbol(DependendSymbol): Line 2472  class Exp_Symbol(DependendSymbol):
2472        @type format: C{str}        @type format: C{str}
2473        @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.
2474        @rtype: C{str}        @rtype: C{str}
2475        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2476        """        """
2477        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2478            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2439  def sqrt(arg): Line 2524  def sqrt(arg):
2524    
2525     @param arg: argument     @param arg: argument
2526     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2527     @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.
2528     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2529     """     """
2530     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2477  class Sqrt_Symbol(DependendSymbol): Line 2562  class Sqrt_Symbol(DependendSymbol):
2562        @type format: C{str}        @type format: C{str}
2563        @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.
2564        @rtype: C{str}        @rtype: C{str}
2565        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2566        """        """
2567        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2568            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2529  def log(arg): Line 2614  def log(arg):
2614    
2615     @param arg: argument     @param arg: argument
2616     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2617     @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.
2618     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2619     """     """
2620     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2567  class Log_Symbol(DependendSymbol): Line 2652  class Log_Symbol(DependendSymbol):
2652        @type format: C{str}        @type format: C{str}
2653        @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.
2654        @rtype: C{str}        @rtype: C{str}
2655        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2656        """        """
2657        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2658            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2619  def sign(arg): Line 2704  def sign(arg):
2704    
2705     @param arg: argument     @param arg: argument
2706     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2707     @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.
2708     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2709     """     """
2710     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2667  class Abs_Symbol(DependendSymbol): Line 2752  class Abs_Symbol(DependendSymbol):
2752        @type format: C{str}        @type format: C{str}
2753        @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.
2754        @rtype: C{str}        @rtype: C{str}
2755        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2756        """        """
2757        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2758            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2719  def minval(arg): Line 2804  def minval(arg):
2804    
2805     @param arg: argument     @param arg: argument
2806     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2807     @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.
2808     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2809     """     """
2810     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2760  class Minval_Symbol(DependendSymbol): Line 2845  class Minval_Symbol(DependendSymbol):
2845        @type format: C{str}        @type format: C{str}
2846        @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.
2847        @rtype: C{str}        @rtype: C{str}
2848        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2849        """        """
2850        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2851            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2796  def maxval(arg): Line 2881  def maxval(arg):
2881    
2882     @param arg: argument     @param arg: argument
2883     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2884     @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.
2885     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2886     """     """
2887     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2837  class Maxval_Symbol(DependendSymbol): Line 2922  class Maxval_Symbol(DependendSymbol):
2922        @type format: C{str}        @type format: C{str}
2923        @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.
2924        @rtype: C{str}        @rtype: C{str}
2925        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2926        """        """
2927        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2928            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2873  def length(arg): Line 2958  def length(arg):
2958    
2959     @param arg: argument     @param arg: argument
2960     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2961     @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.
2962     """     """
2963     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2964    
2965    def trace(arg,axis_offset=0):
2966       """
2967       returns the trace of arg which the sum of arg[k,k] over k.
2968    
2969       @param arg: argument
2970       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2971       @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
2972                      C{axis_offset} and axis_offset+1 must be equal.
2973       @type axis_offset: C{int}
2974       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2975       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2976       """
2977       if isinstance(arg,numarray.NumArray):
2978          sh=arg.shape
2979          if len(sh)<2:
2980            raise ValueError,"rank of argument must be greater than 1"
2981          if axis_offset<0 or axis_offset>len(sh)-2:
2982            raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
2983          s1=1
2984          for i in range(axis_offset): s1*=sh[i]
2985          s2=1
2986          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2987          if not sh[axis_offset] == sh[axis_offset+1]:
2988            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2989          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
2990          out=numarray.zeros([s1,s2],numarray.Float64)
2991          for i1 in range(s1):
2992            for i2 in range(s2):
2993                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
2994          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2995          return out
2996       elif isinstance(arg,escript.Data):
2997          if arg.getRank()<2:
2998            raise ValueError,"rank of argument must be greater than 1"
2999          if axis_offset<0 or axis_offset>arg.getRank()-2:
3000            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3001          s=list(arg.getShape())        
3002          if not s[axis_offset] == s[axis_offset+1]:
3003            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3004          return arg._trace(axis_offset)
3005       elif isinstance(arg,float):
3006          raise TypeError,"illegal argument type float."
3007       elif isinstance(arg,int):
3008          raise TypeError,"illegal argument type int."
3009       elif isinstance(arg,Symbol):
3010          return Trace_Symbol(arg,axis_offset)
3011       else:
3012          raise TypeError,"Unknown argument type."
3013    
3014    class Trace_Symbol(DependendSymbol):
3015       """
3016       L{Symbol} representing the result of the trace function
3017       """
3018       def __init__(self,arg,axis_offset=0):
3019          """
3020          initialization of trace L{Symbol} with argument arg
3021          @param arg: argument of function
3022          @type arg: L{Symbol}.
3023          @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
3024                      C{axis_offset} and axis_offset+1 must be equal.
3025          @type axis_offset: C{int}
3026          """
3027          if arg.getRank()<2:
3028            raise ValueError,"rank of argument must be greater than 1"
3029          if axis_offset<0 or axis_offset>arg.getRank()-2:
3030            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3031          s=list(arg.getShape())        
3032          if not s[axis_offset] == s[axis_offset+1]:
3033            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3034          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3035    
3036       def getMyCode(self,argstrs,format="escript"):
3037          """
3038          returns a program code that can be used to evaluate the symbol.
3039    
3040          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3041          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3042          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3043          @type format: C{str}
3044          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3045          @rtype: C{str}
3046          @raise NotImplementedError: if the requested format is not available
3047          """
3048          if format=="escript" or format=="str"  or format=="text":
3049             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3050          else:
3051             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3052    
3053       def substitute(self,argvals):
3054          """
3055          assigns new values to symbols in the definition of the symbol.
3056          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3057    
3058          @param argvals: new values assigned to symbols
3059          @type argvals: C{dict} with keywords of type L{Symbol}.
3060          @return: result of the substitution process. Operations are executed as much as possible.
3061          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3062          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3063          """
3064          if argvals.has_key(self):
3065             arg=argvals[self]
3066             if self.isAppropriateValue(arg):
3067                return arg
3068             else:
3069                raise TypeError,"%s: new value is not appropriate."%str(self)
3070          else:
3071             arg=self.getSubstitutedArguments(argvals)
3072             return trace(arg[0],axis_offset=arg[1])
3073    
3074       def diff(self,arg):
3075          """
3076          differential of this object
3077    
3078          @param arg: the derivative is calculated with respect to arg
3079          @type arg: L{escript.Symbol}
3080          @return: derivative with respect to C{arg}
3081          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3082          """
3083          if arg==self:
3084             return identity(self.getShape())
3085          else:
3086             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3087    
3088    def transpose(arg,axis_offset=None):
3089       """
3090       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3091    
3092       @param arg: argument
3093       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3094       @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.
3095                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3096       @type axis_offset: C{int}
3097       @return: transpose of arg
3098       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3099       """
3100       if isinstance(arg,numarray.NumArray):
3101          if axis_offset==None: axis_offset=int(arg.rank/2)
3102          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3103       elif isinstance(arg,escript.Data):
3104          r=arg.getRank()
3105          if axis_offset==None: axis_offset=int(r/2)
3106          if axis_offset<0 or axis_offset>r:
3107            raise ValueError,"axis_offset must be between 0 and %s"%r
3108          return arg._transpose(axis_offset)
3109       elif isinstance(arg,float):
3110          if not ( axis_offset==0 or axis_offset==None):
3111            raise ValueError,"axis_offset must be 0 for float argument"
3112          return arg
3113       elif isinstance(arg,int):
3114          if not ( axis_offset==0 or axis_offset==None):
3115            raise ValueError,"axis_offset must be 0 for int argument"
3116          return float(arg)
3117       elif isinstance(arg,Symbol):
3118          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3119          return Transpose_Symbol(arg,axis_offset)
3120       else:
3121          raise TypeError,"Unknown argument type."
3122    
3123    class Transpose_Symbol(DependendSymbol):
3124       """
3125       L{Symbol} representing the result of the transpose function
3126       """
3127       def __init__(self,arg,axis_offset=None):
3128          """
3129          initialization of transpose L{Symbol} with argument arg
3130    
3131          @param arg: argument of function
3132          @type arg: L{Symbol}.
3133          @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.
3134                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3135          @type axis_offset: C{int}
3136          """
3137          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3138          if axis_offset<0 or axis_offset>arg.getRank():
3139            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3140          s=arg.getShape()
3141          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3142    
3143       def getMyCode(self,argstrs,format="escript"):
3144          """
3145          returns a program code that can be used to evaluate the symbol.
3146    
3147          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3148          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3149          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3150          @type format: C{str}
3151          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3152          @rtype: C{str}
3153          @raise NotImplementedError: if the requested format is not available
3154          """
3155          if format=="escript" or format=="str"  or format=="text":
3156             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3157          else:
3158             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3159    
3160       def substitute(self,argvals):
3161          """
3162          assigns new values to symbols in the definition of the symbol.
3163          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3164    
3165          @param argvals: new values assigned to symbols
3166          @type argvals: C{dict} with keywords of type L{Symbol}.
3167          @return: result of the substitution process. Operations are executed as much as possible.
3168          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3169          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3170          """
3171          if argvals.has_key(self):
3172             arg=argvals[self]
3173             if self.isAppropriateValue(arg):
3174                return arg
3175             else:
3176                raise TypeError,"%s: new value is not appropriate."%str(self)
3177          else:
3178             arg=self.getSubstitutedArguments(argvals)
3179             return transpose(arg[0],axis_offset=arg[1])
3180    
3181       def diff(self,arg):
3182          """
3183          differential of this object
3184    
3185          @param arg: the derivative is calculated with respect to arg
3186          @type arg: L{escript.Symbol}
3187          @return: derivative with respect to C{arg}
3188          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3189          """
3190          if arg==self:
3191             return identity(self.getShape())
3192          else:
3193             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3194    
3195    def swap_axes(arg,axis0=0,axis1=1):
3196       """
3197       returns the swap of arg by swaping the components axis0 and axis1
3198    
3199       @param arg: argument
3200       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3201       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3202       @type axis0: C{int}
3203       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3204       @type axis1: C{int}
3205       @return: C{arg} with swaped components
3206       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3207       """
3208       if axis0 > axis1:
3209          axis0,axis1=axis1,axis0
3210       if isinstance(arg,numarray.NumArray):
3211          return numarray.swapaxes(arg,axis0,axis1)
3212       elif isinstance(arg,escript.Data):
3213          return arg._swap_axes(axis0,axis1)
3214       elif isinstance(arg,float):
3215          raise TyepError,"float argument is not supported."
3216       elif isinstance(arg,int):
3217          raise TyepError,"int argument is not supported."
3218       elif isinstance(arg,Symbol):
3219          return SwapAxes_Symbol(arg,axis0,axis1)
3220       else:
3221          raise TypeError,"Unknown argument type."
3222    
3223    class SwapAxes_Symbol(DependendSymbol):
3224       """
3225       L{Symbol} representing the result of the swap function
3226       """
3227       def __init__(self,arg,axis0=0,axis1=1):
3228          """
3229          initialization of swap L{Symbol} with argument arg
3230    
3231          @param arg: argument
3232          @type arg: L{Symbol}.
3233          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3234          @type axis0: C{int}
3235          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3236          @type axis1: C{int}
3237          """
3238          if arg.getRank()<2:
3239             raise ValueError,"argument must have at least rank 2."
3240          if axis0<0 or axis0>arg.getRank()-1:
3241             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3242          if axis1<0 or axis1>arg.getRank()-1:
3243             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3244          if axis0 == axis1:
3245             raise ValueError,"axis indices must be different."
3246          if axis0 > axis1:
3247             axis0,axis1=axis1,axis0
3248          s=arg.getShape()
3249          s_out=[]
3250          for i in range(len(s)):
3251             if i == axis0:
3252                s_out.append(s[axis1])
3253             elif i == axis1:
3254                s_out.append(s[axis0])
3255             else:
3256                s_out.append(s[i])
3257          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3258    
3259       def getMyCode(self,argstrs,format="escript"):
3260          """
3261          returns a program code that can be used to evaluate the symbol.
3262    
3263          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3264          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3265          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3266          @type format: C{str}
3267          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3268          @rtype: C{str}
3269          @raise NotImplementedError: if the requested format is not available
3270          """
3271          if format=="escript" or format=="str"  or format=="text":
3272             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3273          else:
3274             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3275    
3276       def substitute(self,argvals):
3277          """
3278          assigns new values to symbols in the definition of the symbol.
3279          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3280    
3281          @param argvals: new values assigned to symbols
3282          @type argvals: C{dict} with keywords of type L{Symbol}.
3283          @return: result of the substitution process. Operations are executed as much as possible.
3284          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3285          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3286          """
3287          if argvals.has_key(self):
3288             arg=argvals[self]
3289             if self.isAppropriateValue(arg):
3290                return arg
3291             else:
3292                raise TypeError,"%s: new value is not appropriate."%str(self)
3293          else:
3294             arg=self.getSubstitutedArguments(argvals)
3295             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3296    
3297       def diff(self,arg):
3298          """
3299          differential of this object
3300    
3301          @param arg: the derivative is calculated with respect to arg
3302          @type arg: L{escript.Symbol}
3303          @return: derivative with respect to C{arg}
3304          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3305          """
3306          if arg==self:
3307             return identity(self.getShape())
3308          else:
3309             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3310    
3311    def symmetric(arg):
3312        """
3313        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3314    
3315        @param arg: square matrix. Must have rank 2 or 4 and be square.
3316        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3317        @return: symmetric part of arg
3318        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3319        """
3320        if isinstance(arg,numarray.NumArray):
3321          if arg.rank==2:
3322            if not (arg.shape[0]==arg.shape[1]):
3323               raise ValueError,"argument must be square."
3324          elif arg.rank==4:
3325            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3326               raise ValueError,"argument must be square."
3327          else:
3328            raise ValueError,"rank 2 or 4 is required."
3329          return (arg+transpose(arg))/2
3330        elif isinstance(arg,escript.Data):
3331          if arg.getRank()==2:
3332            if not (arg.getShape()[0]==arg.getShape()[1]):
3333               raise ValueError,"argument must be square."
3334            return arg._symmetric()
3335          elif arg.getRank()==4:
3336            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3337               raise ValueError,"argument must be square."
3338            return arg._symmetric()
3339          else:
3340            raise ValueError,"rank 2 or 4 is required."
3341        elif isinstance(arg,float):
3342          return arg
3343        elif isinstance(arg,int):
3344          return float(arg)
3345        elif isinstance(arg,Symbol):
3346          if arg.getRank()==2:
3347            if not (arg.getShape()[0]==arg.getShape()[1]):
3348               raise ValueError,"argument must be square."
3349          elif arg.getRank()==4:
3350            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3351               raise ValueError,"argument must be square."
3352          else:
3353            raise ValueError,"rank 2 or 4 is required."
3354          return (arg+transpose(arg))/2
3355        else:
3356          raise TypeError,"symmetric: Unknown argument type."
3357    
3358    def nonsymmetric(arg):
3359        """
3360        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3361    
3362        @param arg: square matrix. Must have rank 2 or 4 and be square.
3363        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3364        @return: nonsymmetric part of arg
3365        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3366        """
3367        if isinstance(arg,numarray.NumArray):
3368          if arg.rank==2:
3369            if not (arg.shape[0]==arg.shape[1]):
3370               raise ValueError,"nonsymmetric: argument must be square."
3371          elif arg.rank==4:
3372            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3373               raise ValueError,"nonsymmetric: argument must be square."
3374          else:
3375            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3376          return (arg-transpose(arg))/2
3377        elif isinstance(arg,escript.Data):
3378          if arg.getRank()==2:
3379            if not (arg.getShape()[0]==arg.getShape()[1]):
3380               raise ValueError,"argument must be square."
3381            return arg._nonsymmetric()
3382          elif arg.getRank()==4:
3383            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3384               raise ValueError,"argument must be square."
3385            return arg._nonsymmetric()
3386          else:
3387            raise ValueError,"rank 2 or 4 is required."
3388        elif isinstance(arg,float):
3389          return arg
3390        elif isinstance(arg,int):
3391          return float(arg)
3392        elif isinstance(arg,Symbol):
3393          if arg.getRank()==2:
3394            if not (arg.getShape()[0]==arg.getShape()[1]):
3395               raise ValueError,"nonsymmetric: argument must be square."
3396          elif arg.getRank()==4:
3397            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3398               raise ValueError,"nonsymmetric: argument must be square."
3399          else:
3400            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3401          return (arg-transpose(arg))/2
3402        else:
3403          raise TypeError,"nonsymmetric: Unknown argument type."
3404    
3405    def inverse(arg):
3406        """
3407        returns the inverse of the square matrix arg.
3408    
3409        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3410        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3411        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3412        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3413        @note: for L{escript.Data} objects the dimension is restricted to 3.
3414        """
3415        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3416        if isinstance(arg,numarray.NumArray):
3417          return numarray.linear_algebra.inverse(arg)
3418        elif isinstance(arg,escript.Data):
3419          return escript_inverse(arg)
3420        elif isinstance(arg,float):
3421          return 1./arg
3422        elif isinstance(arg,int):
3423          return 1./float(arg)
3424        elif isinstance(arg,Symbol):
3425          return Inverse_Symbol(arg)
3426        else:
3427          raise TypeError,"inverse: Unknown argument type."
3428    
3429    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3430          "arg is a Data objects!!!"
3431          if not arg.getRank()==2:
3432            raise ValueError,"escript_inverse: argument must have rank 2"
3433          s=arg.getShape()      
3434          if not s[0] == s[1]:
3435            raise ValueError,"escript_inverse: argument must be a square matrix."
3436          out=escript.Data(0.,s,arg.getFunctionSpace())
3437          if s[0]==1:
3438              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3439                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3440              out[0,0]=1./arg[0,0]
3441          elif s[0]==2:
3442              A11=arg[0,0]
3443              A12=arg[0,1]
3444              A21=arg[1,0]
3445              A22=arg[1,1]
3446              D = A11*A22-A12*A21
3447              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3448                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3449              D=1./D
3450              out[0,0]= A22*D
3451              out[1,0]=-A21*D
3452              out[0,1]=-A12*D
3453              out[1,1]= A11*D
3454          elif s[0]==3:
3455              A11=arg[0,0]
3456              A21=arg[1,0]
3457              A31=arg[2,0]
3458              A12=arg[0,1]
3459              A22=arg[1,1]
3460              A32=arg[2,1]
3461              A13=arg[0,2]
3462              A23=arg[1,2]
3463              A33=arg[2,2]
3464              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3465              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3466                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3467              D=1./D
3468              out[0,0]=(A22*A33-A23*A32)*D
3469              out[1,0]=(A31*A23-A21*A33)*D
3470              out[2,0]=(A21*A32-A31*A22)*D
3471              out[0,1]=(A13*A32-A12*A33)*D
3472              out[1,1]=(A11*A33-A31*A13)*D
3473              out[2,1]=(A12*A31-A11*A32)*D
3474              out[0,2]=(A12*A23-A13*A22)*D
3475              out[1,2]=(A13*A21-A11*A23)*D
3476              out[2,2]=(A11*A22-A12*A21)*D
3477          else:
3478             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3479          return out
3480    
3481    class Inverse_Symbol(DependendSymbol):
3482       """
3483       L{Symbol} representing the result of the inverse function
3484       """
3485       def __init__(self,arg):
3486          """
3487          initialization of inverse L{Symbol} with argument arg
3488          @param arg: argument of function
3489          @type arg: L{Symbol}.
3490          """
3491          if not arg.getRank()==2:
3492            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3493          s=arg.getShape()
3494          if not s[0] == s[1]:
3495            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3496          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3497    
3498       def getMyCode(self,argstrs,format="escript"):
3499          """
3500          returns a program code that can be used to evaluate the symbol.
3501    
3502          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3503          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3504          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3505          @type format: C{str}
3506          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3507          @rtype: C{str}
3508          @raise NotImplementedError: if the requested format is not available
3509          """
3510          if format=="escript" or format=="str"  or format=="text":
3511             return "inverse(%s)"%argstrs[0]
3512          else:
3513             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3514    
3515       def substitute(self,argvals):
3516          """
3517          assigns new values to symbols in the definition of the symbol.
3518          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3519    
3520          @param argvals: new values assigned to symbols
3521          @type argvals: C{dict} with keywords of type L{Symbol}.
3522          @return: result of the substitution process. Operations are executed as much as possible.
3523          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3524          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3525          """
3526          if argvals.has_key(self):
3527             arg=argvals[self]
3528             if self.isAppropriateValue(arg):
3529                return arg
3530             else:
3531                raise TypeError,"%s: new value is not appropriate."%str(self)
3532          else:
3533             arg=self.getSubstitutedArguments(argvals)
3534             return inverse(arg[0])
3535    
3536       def diff(self,arg):
3537          """
3538          differential of this object
3539    
3540          @param arg: the derivative is calculated with respect to arg
3541          @type arg: L{escript.Symbol}
3542          @return: derivative with respect to C{arg}
3543          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3544          """
3545          if arg==self:
3546             return identity(self.getShape())
3547          else:
3548             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3549    
3550    def eigenvalues(arg):
3551        """
3552        returns the eigenvalues of the square matrix arg.
3553    
3554        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3555                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3556        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3557        @return: the eigenvalues in increasing order.
3558        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3559        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3560        """
3561        if isinstance(arg,numarray.NumArray):
3562          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3563          out.sort()
3564          return out
3565        elif isinstance(arg,escript.Data):
3566          return arg._eigenvalues()
3567        elif isinstance(arg,Symbol):
3568          if not arg.getRank()==2:
3569            raise ValueError,"eigenvalues: argument must have rank 2"
3570          s=arg.getShape()      
3571          if not s[0] == s[1]:
3572            raise ValueError,"eigenvalues: argument must be a square matrix."
3573          if s[0]==1:
3574              return arg[0]
3575          elif s[0]==2:
3576              arg1=symmetric(arg)
3577              A11=arg1[0,0]
3578              A12=arg1[0,1]
3579              A22=arg1[1,1]
3580              trA=(A11+A22)/2.
3581              A11-=trA
3582              A22-=trA
3583              s=sqrt(A12**2-A11*A22)
3584              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3585          elif s[0]==3:
3586              arg1=symmetric(arg)
3587              A11=arg1[0,0]
3588              A12=arg1[0,1]
3589              A22=arg1[1,1]
3590              A13=arg1[0,2]
3591              A23=arg1[1,2]
3592              A33=arg1[2,2]
3593              trA=(A11+A22+A33)/3.
3594              A11-=trA
3595              A22-=trA
3596              A33-=trA
3597              A13_2=A13**2
3598              A23_2=A23**2
3599              A12_2=A12**2
3600              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3601              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3602              sq_p=sqrt(p/3.)
3603              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
3604              sq_p*=2.
3605              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3606               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3607               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3608              return trA+sq_p*f
3609          else:
3610             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3611        elif isinstance(arg,float):
3612          return arg
3613        elif isinstance(arg,int):
3614          return float(arg)
3615        else:
3616          raise TypeError,"eigenvalues: Unknown argument type."
3617    
3618    def eigenvalues_and_eigenvectors(arg):
3619        """
3620        returns the eigenvalues and eigenvectors of the square matrix arg.
3621    
3622        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3623                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3624        @type arg: L{escript.Data}
3625        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3626                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3627                 the eigenvector coresponding to the i-th eigenvalue.
3628        @rtype: L{tuple} of L{escript.Data}.
3629        @note: The dimension is restricted to 3.
3630        """
3631        if isinstance(arg,numarray.NumArray):
3632          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3633        elif isinstance(arg,escript.Data):
3634          return arg._eigenvalues_and_eigenvectors()
3635        elif isinstance(arg,Symbol):
3636          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3637        elif isinstance(arg,float):
3638          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3639        elif isinstance(arg,int):
3640          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3641        else:
3642          raise TypeError,"eigenvalues: Unknown argument type."
3643  #=======================================================  #=======================================================
3644  #  Binary operations:  #  Binary operations:
3645  #=======================================================  #=======================================================
# Line 2920  class Add_Symbol(DependendSymbol): Line 3683  class Add_Symbol(DependendSymbol):
3683         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3684         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3685         """         """
3686         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3687         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3688         if not sh0==sh1:         if not sh0==sh1:
3689            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3690         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 2936  class Add_Symbol(DependendSymbol): Line 3699  class Add_Symbol(DependendSymbol):
3699        @type format: C{str}        @type format: C{str}
3700        @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.
3701        @rtype: C{str}        @rtype: C{str}
3702        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3703        """        """
3704        if format=="str" or format=="text":        if format=="str" or format=="text":
3705           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 2995  def mult(arg0,arg1): Line 3758  def mult(arg0,arg1):
3758         """         """
3759         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3760         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3761            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3762         else:         else:
3763            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3764                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3019  class Mult_Symbol(DependendSymbol): Line 3782  class Mult_Symbol(DependendSymbol):
3782         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3783         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3784         """         """
3785         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3786         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3787         if not sh0==sh1:         if not sh0==sh1:
3788            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3789         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3035  class Mult_Symbol(DependendSymbol): Line 3798  class Mult_Symbol(DependendSymbol):
3798        @type format: C{str}        @type format: C{str}
3799        @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.
3800        @rtype: C{str}        @rtype: C{str}
3801        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3802        """        """
3803        if format=="str" or format=="text":        if format=="str" or format=="text":
3804           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3095  def quotient(arg0,arg1): Line 3858  def quotient(arg0,arg1):
3858         """         """
3859         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3860         if testForZero(args[0]):         if testForZero(args[0]):
3861            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3862         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3863            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3864               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3124  class Quotient_Symbol(DependendSymbol): Line 3887  class Quotient_Symbol(DependendSymbol):
3887         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3888         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3889         """         """
3890         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3891         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3892         if not sh0==sh1:         if not sh0==sh1:
3893            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3894         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3140  class Quotient_Symbol(DependendSymbol): Line 3903  class Quotient_Symbol(DependendSymbol):
3903        @type format: C{str}        @type format: C{str}
3904        @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.
3905        @rtype: C{str}        @rtype: C{str}
3906        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3907        """        """
3908        if format=="str" or format=="text":        if format=="str" or format=="text":
3909           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3201  def power(arg0,arg1): Line 3964  def power(arg0,arg1):
3964         """         """
3965         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3966         if testForZero(args[0]):         if testForZero(args[0]):
3967            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3968         elif testForZero(args[1]):         elif testForZero(args[1]):
3969            return numarray.ones(args[0],numarray.Float)            return numarray.ones(getShape(args[1]),numarray.Float64)
3970         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3971            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3972         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):
# Line 3226  class Power_Symbol(DependendSymbol): Line 3989  class Power_Symbol(DependendSymbol):
3989         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3990         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3991         """         """
3992         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3993         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3994         if not sh0==sh1:         if not sh0==sh1:
3995            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
3996         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3244  class Power_Symbol(DependendSymbol): Line 4007  class Power_Symbol(DependendSymbol):
4007        @type format: C{str}        @type format: C{str}
4008        @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.
4009        @rtype: C{str}        @rtype: C{str}
4010        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4011        """        """
4012        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4013           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3304  def maximum(*args): Line 4067  def maximum(*args):
4067         if out==None:         if out==None:
4068            out=a            out=a
4069         else:         else:
4070            m=whereNegative(out-a)            diff=add(a,-out)
4071            out=m*a+(1.-m)*out            out=add(out,mult(wherePositive(diff),diff))
4072      return out      return out
4073        
4074  def minimum(*arg):  def minimum(*args):
4075      """      """
4076      the minimum over arguments args      the minimum over arguments args
4077    
# Line 3322  def minimum(*arg): Line 4085  def minimum(*arg):
4085         if out==None:         if out==None:
4086            out=a            out=a
4087         else:         else:
4088            m=whereNegative(out-a)            diff=add(a,-out)
4089            out=m*out+(1.-m)*a            out=add(out,mult(whereNegative(diff),diff))
4090      return out      return out
4091    
4092    def clip(arg,minval=None,maxval=None):
4093        """
4094        cuts the values of arg between minval and maxval
4095    
4096        @param arg: argument
4097        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4098        @param minval: lower range. If None no lower range is applied
4099        @type minval: C{float} or C{None}
4100        @param maxval: upper range. If None no upper range is applied
4101        @type maxval: C{float} or C{None}
4102        @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4103                 less then maxval are unchanged.
4104        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4105        @raise ValueError: if minval>maxval
4106        """
4107        if not minval==None and not maxval==None:
4108           if minval>maxval:
4109              raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4110        if minval == None:
4111            tmp=arg
4112        else:
4113            tmp=maximum(minval,arg)
4114        if maxval == None:
4115            return tmp
4116        else:
4117            return minimum(tmp,maxval)
4118    
4119        
4120  def inner(arg0,arg1):  def inner(arg0,arg1):
4121      """      """
# Line 3340  def inner(arg0,arg1): Line 4131  def inner(arg0,arg1):
4131      @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}
4132      @param arg1: second argument      @param arg1: second argument
4133      @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}
4134      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4135      @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
4136      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4137      """      """
4138      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4139      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4140      if not sh0==sh1:      if not sh0==sh1:
4141          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4142      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4143    
4144    def outer(arg0,arg1):
4145        """
4146        the outer product of the two argument:
4147    
4148        out[t,s]=arg0[t]*arg1[s]
4149    
4150        where
4151    
4152            - s runs through arg0.Shape
4153            - t runs through arg1.Shape
4154    
4155        @param arg0: first argument
4156        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4157        @param arg1: second argument
4158        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4159        @return: the outer product of arg0 and arg1 at each data point
4160        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4161        """
4162        return generalTensorProduct(arg0,arg1,axis_offset=0)
4163    
4164  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4165      """      """
4166        see L{matrix_mult}
4167        """
4168        return matrix_mult(arg0,arg1)
4169    
4170    def matrix_mult(arg0,arg1):
4171        """
4172      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4173    
4174      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4175    
4176            or      or
4177    
4178      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4179    
4180      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.
4181    
4182      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4183      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3370  def matrixmult(arg0,arg1): Line 4187  def matrixmult(arg0,arg1):
4187      @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
4188      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4189      """      """
4190      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4191      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4192      if not len(sh0)==2 :      if not len(sh0)==2 :
4193          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4194      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4195          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4196      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4197    
4198  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4199      """      """
4200      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  
4201      """      """
4202      return generalTensorProduct(arg0,arg1,offset=0)      return tensor_mult(arg0,arg1)
4203    
4204    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4205      """      """
4206      the tensor product of the two argument:      the tensor product of the two argument:
   
4207            
4208      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4209    
4210      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4211    
4212                   or      or
4213    
4214      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4215    
# Line 3415  def tensormult(arg0,arg1): Line 4218  def tensormult(arg0,arg1):
4218    
4219      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]
4220                                
4221                   or      or
4222    
4223      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]
4224    
4225                   or      or
4226    
4227      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]
4228    
4229      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  
4230      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.
4231    
4232      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4233      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3433  def tensormult(arg0,arg1): Line 4236  def tensormult(arg0,arg1):
4236      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4237      @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
4238      """      """
4239      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4240      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4241      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4242         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4243      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):
4244         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4245      else:      else:
4246          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4247    
4248  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4249      """      """
4250      generalized tensor product      generalized tensor product
4251    
4252      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4253    
4254      where s runs through arg0.Shape[:arg0.Rank-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[offset:]  
4255    
4256      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]
4257      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4258            - t runs through arg1.Shape[axis_offset:]
4259    
4260      @param arg0: first argument      @param arg0: first argument
4261      @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}
4262      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4263      @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}
4264      @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.
4265      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
# Line 3467  def generalTensorProduct(arg0,arg1,offse Line 4269  def generalTensorProduct(arg0,arg1,offse
4269      # 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
4270      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4271         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4272             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4273         else:         else:
4274             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4275                 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)
4276             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4277             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4278             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4279             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4280             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4281             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4282             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4283             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4284             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4285             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4286             for i0 in range(d0):             for i0 in range(d0):
4287                      for i1 in range(d1):                      for i1 in range(d1):
4288                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4289             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4290             return out             return out
4291      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4292         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4293             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4294         else:         else:
4295             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)
4296      else:            else:      
4297         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4298                                    
4299  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4300     """     """
4301     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4302     """     """
4303     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4304         """         """
4305         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4306    
4307         @param arg0: numerator         @param arg0: first argument
4308         @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}.
4309         @param arg1: denominator         @param arg1: second argument
4310         @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}.
4311         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4312         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4313         """         """
4314         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4315         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4316         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4317         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4318         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4319         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4320         if not sh01==sh10:         if not sh01==sh10:
4321             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)
4322         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])
4323    
4324     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4325        """        """
# Line 3529  class GeneralTensorProduct_Symbol(Depend Line 4331  class GeneralTensorProduct_Symbol(Depend
4331        @type format: C{str}        @type format: C{str}
4332        @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.
4333        @rtype: C{str}        @rtype: C{str}
4334        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4335        """        """
4336        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4337           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])
4338        else:        else:
4339           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)
4340    
# Line 3557  class GeneralTensorProduct_Symbol(Depend Line 4359  class GeneralTensorProduct_Symbol(Depend
4359           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4360           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4361    
4362  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4363      "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!!!"
4364      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4365      shape0=arg0.getShape()[:arg0.getRank()-offset]  
4366      shape01=arg0.getShape()[arg0.getRank()-offset:]  def transposed_matrix_mult(arg0,arg1):
4367      shape10=arg1.getShape()[:offset]      """
4368      shape1=arg1.getShape()[offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4369      if not shape01==shape10:  
4370          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]
4371    
4372      # create return value:      or
4373      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())  
4374      #      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4375      s0=[[]]  
4376      for k in shape0:      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4377            s=[]  
4378            for j in s0:      The first dimension of arg0 and arg1 must match.
4379                  for i in range(k): s.append(j+[slice(i,i)])  
4380            s0=s      @param arg0: first argument of rank 2
4381      s1=[[]]      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4382      for k in shape1:      @param arg1: second argument of at least rank 1
4383            s=[]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4384            for j in s1:      @return: the product of the transposed of arg0 and arg1 at each data point
4385                  for i in range(k): s.append(j+[slice(i,i)])      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4386            s1=s      @raise ValueError: if the shapes of the arguments are not appropriate
4387      s01=[[]]      """
4388      for k in shape01:      sh0=getShape(arg0)
4389            s=[]      sh1=getShape(arg1)
4390            for j in s01:      if not len(sh0)==2 :
4391                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"first argument must have rank 2"
4392            s01=s      if not len(sh1)==2 and not len(sh1)==1:
4393            raise ValueError,"second argument must have rank 1 or 2"
4394      for i0 in s0:      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4395         for i1 in s1:  
4396           s=escript.Scalar(0.,arg0.getFunctionSpace())  def transposed_tensor_mult(arg0,arg1):
4397           for i01 in s01:      """
4398              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      the tensor product of the transposed of the first and the second argument
4399           out.__setitem__(tuple(i0+i1),s)      
4400      return out      for arg0 of rank 2 this is
4401    
4402        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4403    
4404        or
4405    
4406        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4407    
4408      
4409        and for arg0 of rank 4 this is
4410    
4411        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4412                  
4413        or
4414    
4415        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4416    
4417        or
4418    
4419        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4420    
4421        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4422        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4423    
4424        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4425    
4426        @param arg0: first argument of rank 2 or 4
4427        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4428        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4429        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4430        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4431        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4432        """
4433        sh0=getShape(arg0)
4434        sh1=getShape(arg1)
4435        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4436           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4437        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4438           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4439        else:
4440            raise ValueError,"first argument must have rank 2 or 4"
4441    
4442    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4443        """
4444        generalized tensor product of transposed of arg0 and arg1:
4445    
4446        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4447    
4448        where
4449    
4450            - s runs through arg0.Shape[axis_offset:]
4451            - r runs trough arg0.Shape[:axis_offset]
4452            - t runs through arg1.Shape[axis_offset:]
4453    
4454        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4455        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4456    
4457        @param arg0: first argument
4458        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4459        @param arg1: second argument
4460        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4461        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4462        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4463        """
4464        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4465        arg0,arg1=matchType(arg0,arg1)
4466        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4467        if isinstance(arg0,numarray.NumArray):
4468           if isinstance(arg1,Symbol):
4469               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4470           else:
4471               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4472                   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)
4473               arg0_c=arg0.copy()
4474               arg1_c=arg1.copy()
4475               sh0,sh1=arg0.shape,arg1.shape
4476               d0,d1,d01=1,1,1
4477               for i in sh0[axis_offset:]: d0*=i
4478               for i in sh1[axis_offset:]: d1*=i
4479               for i in sh0[:axis_offset]: d01*=i
4480               arg0_c.resize((d01,d0))
4481               arg1_c.resize((d01,d1))
4482               out=numarray.zeros((d0,d1),numarray.Float64)
4483               for i0 in range(d0):
4484                        for i1 in range(d1):
4485                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4486               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4487               return out
4488        elif isinstance(arg0,escript.Data):
4489           if isinstance(arg1,Symbol):
4490               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4491           else:
4492               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4493        else:      
4494           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4495                    
4496    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4497       """
4498       Symbol representing the general tensor product of the transposed of arg0 and arg1
4499       """
4500       def __init__(self,arg0,arg1,axis_offset=0):
4501           """
4502           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4503    
4504           @param arg0: first argument
4505           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4506           @param arg1: second argument
4507           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4508           @raise ValueError: inconsistent dimensions of arguments.
4509           @note: if both arguments have a spatial dimension, they must equal.
4510           """
4511           sh_arg0=getShape(arg0)
4512           sh_arg1=getShape(arg1)
4513           sh01=sh_arg0[:axis_offset]
4514           sh10=sh_arg1[:axis_offset]
4515           sh0=sh_arg0[axis_offset:]
4516           sh1=sh_arg1[axis_offset:]
4517           if not sh01==sh10:
4518               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)
4519           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4520    
4521       def getMyCode(self,argstrs,format="escript"):
4522          """
4523          returns a program code that can be used to evaluate the symbol.
4524    
4525          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4526          @type argstrs: C{list} of length 2 of C{str}.
4527          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4528          @type format: C{str}
4529          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4530          @rtype: C{str}
4531          @raise NotImplementedError: if the requested format is not available
4532          """
4533          if format=="escript" or format=="str" or format=="text":
4534             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4535          else:
4536             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4537    
4538       def substitute(self,argvals):
4539          """
4540          assigns new values to symbols in the definition of the symbol.
4541          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4542    
4543          @param argvals: new values assigned to symbols
4544          @type argvals: C{dict} with keywords of type L{Symbol}.
4545          @return: result of the substitution process. Operations are executed as much as possible.
4546          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4547          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4548          """
4549          if argvals.has_key(self):
4550             arg=argvals[self]
4551             if self.isAppropriateValue(arg):
4552                return arg
4553             else:
4554                raise TypeError,"%s: new value is not appropriate."%str(self)
4555          else:
4556             args=self.getSubstitutedArguments(argvals)
4557             return generalTransposedTensorProduct(args[0],args[1],args[2])
4558    
4559    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4560        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4561        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4562    
4563    def matrix_transposed_mult(arg0,arg1):
4564        """
4565        matrix-transposed(matrix) product of the two argument:
4566    
4567        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4568    
4569        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4570    
4571        The last dimensions of arg0 and arg1 must match.
4572    
4573        @param arg0: first argument of rank 2
4574        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4575        @param arg1: second argument of rank 2
4576        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4577        @return: the product of arg0 and the transposed of arg1 at each data point
4578        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4579        @raise ValueError: if the shapes of the arguments are not appropriate
4580        """
4581        sh0=getShape(arg0)
4582        sh1=getShape(arg1)
4583        if not len(sh0)==2 :
4584            raise ValueError,"first argument must have rank 2"
4585        if not len(sh1)==2 and not len(sh1)==1:
4586            raise ValueError,"second argument must have rank 1 or 2"
4587        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4588    
4589    def tensor_transposed_mult(arg0,arg1):
4590        """
4591        the tensor product of the first and the transpose of the second argument
4592        
4593        for arg0 of rank 2 this is
4594    
4595        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4596    
4597        and for arg0 of rank 4 this is
4598    
4599        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4600                  
4601        or
4602    
4603        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4604    
4605        In the first case the the second dimension of arg0 and arg1 must match and  
4606        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4607    
4608        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4609    
4610        @param arg0: first argument of rank 2 or 4
4611        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4612        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4613        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4614        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4615        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4616        """
4617        sh0=getShape(arg0)
4618        sh1=getShape(arg1)
4619        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4620           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4621        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4622           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4623        else:
4624            raise ValueError,"first argument must have rank 2 or 4"
4625    
4626    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4627        """
4628        generalized tensor product of transposed of arg0 and arg1:
4629    
4630        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4631    
4632        where
4633    
4634            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4635            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4636            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4637    
4638        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4639        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4640    
4641        @param arg0: first argument
4642        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4643        @param arg1: second argument
4644        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4645        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4646        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4647        """
4648        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4649        arg0,arg1=matchType(arg0,arg1)
4650        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4651        if isinstance(arg0,numarray.NumArray):
4652           if isinstance(arg1,Symbol):
4653               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4654           else:
4655               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4656                   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)
4657               arg0_c=arg0.copy()
4658               arg1_c=arg1.copy()
4659               sh0,sh1=arg0.shape,arg1.shape
4660               d0,d1,d01=1,1,1
4661               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4662               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4663               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4664               arg0_c.resize((d0,d01))
4665               arg1_c.resize((d1,d01))
4666               out=numarray.zeros((d0,d1),numarray.Float64)
4667               for i0 in range(d0):
4668                        for i1 in range(d1):
4669                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4670               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4671               return out
4672        elif isinstance(arg0,escript.Data):
4673           if isinstance(arg1,Symbol):
4674               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4675           else:
4676               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4677        else:      
4678           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4679                    
4680    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4681       """
4682       Symbol representing the general tensor product of arg0 and the transpose of arg1
4683       """
4684       def __init__(self,arg0,arg1,axis_offset=0):
4685           """
4686           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4687    
4688           @param arg0: first argument
4689           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4690           @param arg1: second argument
4691           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4692           @raise ValueError: inconsistent dimensions of arguments.
4693           @note: if both arguments have a spatial dimension, they must equal.
4694           """
4695           sh_arg0=getShape(arg0)
4696           sh_arg1=getShape(arg1)
4697           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4698           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4699           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4700           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4701           if not sh01==sh10:
4702               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)
4703           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4704    
4705       def getMyCode(self,argstrs,format="escript"):
4706          """
4707          returns a program code that can be used to evaluate the symbol.
4708    
4709          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4710          @type argstrs: C{list} of length 2 of C{str}.
4711          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4712          @type format: C{str}
4713          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4714          @rtype: C{str}
4715          @raise NotImplementedError: if the requested format is not available
4716          """
4717          if format=="escript" or format=="str" or format=="text":
4718             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4719          else:
4720             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4721    
4722       def substitute(self,argvals):
4723          """
4724          assigns new values to symbols in the definition of the symbol.
4725          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4726    
4727          @param argvals: new values assigned to symbols
4728          @type argvals: C{dict} with keywords of type L{Symbol}.
4729          @return: result of the substitution process. Operations are executed as much as possible.
4730          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4731          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4732          """
4733          if argvals.has_key(self):
4734             arg=argvals[self]
4735             if self.isAppropriateValue(arg):
4736                return arg
4737             else:
4738                raise TypeError,"%s: new value is not appropriate."%str(self)
4739          else:
4740             args=self.getSubstitutedArguments(argvals)
4741             return generalTensorTransposedProduct(args[0],args[1],args[2])
4742    
4743    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4744        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4745        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4746    
4747  #=========================================================  #=========================================================
4748  #   some little helpers  #  functions dealing with spatial dependency
4749  #=========================================================  #=========================================================
4750  def grad(arg,where=None):  def grad(arg,where=None):
4751      """      """
4752      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4753    
4754        If C{g} is the returned object, then
4755    
4756      @param arg:   Data object representing the function which gradient        - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4757                    to be calculated.        - if C{arg} is rank 1 C{g[i,s]} is the derivative of C{arg[i]} with respect to the C{s}-th spatial dimension.
4758          - 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.
4759          - 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.
4760    
4761        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4762        @type arg: L{escript.Data} or L{Symbol}
4763      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4764                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4765        @type where: C{None} or L{escript.FunctionSpace}
4766        @return: gradient of arg.
4767        @rtype: L{escript.Data} or L{Symbol}
4768      """      """
4769      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4770         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3617  def grad(arg,where=None): Line 4774  def grad(arg,where=None):
4774         else:         else:
4775            return arg._grad(where)            return arg._grad(where)
4776      else:      else:
4777        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4778    
4779    class Grad_Symbol(DependendSymbol):
4780       """
4781       L{Symbol} representing the result of the gradient operator
4782       """
4783       def __init__(self,arg,where=None):
4784          """
4785          initialization of gradient L{Symbol} with argument arg
4786          @param arg: argument of function
4787          @type arg: L{Symbol}.
4788          @param where: FunctionSpace in which the gradient will be calculated.
4789                      If not present or C{None} an appropriate default is used.
4790          @type where: C{None} or L{escript.FunctionSpace}
4791          """
4792          d=arg.getDim()
4793          if d==None:
4794             raise ValueError,"argument must have a spatial dimension"
4795          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4796    
4797       def getMyCode(self,argstrs,format="escript"):
4798          """
4799          returns a program code that can be used to evaluate the symbol.
4800    
4801          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4802          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4803          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4804          @type format: C{str}
4805          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4806          @rtype: C{str}
4807          @raise NotImplementedError: if the requested format is not available
4808          """
4809          if format=="escript" or format=="str"  or format=="text":
4810             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4811          else:
4812             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4813    
4814       def substitute(self,argvals):
4815          """
4816          assigns new values to symbols in the definition of the symbol.
4817          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4818    
4819          @param argvals: new values assigned to symbols
4820          @type argvals: C{dict} with keywords of type L{Symbol}.
4821          @return: result of the substitution process. Operations are executed as much as possible.
4822          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4823          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4824          """
4825          if argvals.has_key(self):
4826             arg=argvals[self]
4827             if self.isAppropriateValue(arg):
4828                return arg
4829             else:
4830                raise TypeError,"%s: new value is not appropriate."%str(self)
4831          else:
4832             arg=self.getSubstitutedArguments(argvals)
4833             return grad(arg[0],where=arg[1])
4834    
4835       def diff(self,arg):
4836          """
4837          differential of this object
4838    
4839          @param arg: the derivative is calculated with respect to arg
4840          @type arg: L{escript.Symbol}
4841          @return: derivative with respect to C{arg}
4842          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4843          """
4844          if arg==self:
4845             return identity(self.getShape())
4846          else:
4847             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4848    
4849  def integrate(arg,where=None):  def integrate(arg,where=None):
4850      """      """
4851      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}
4852      its domain.      before integration.
4853    
4854      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4855        @type arg: L{escript.Data} or L{Symbol}
4856      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4857                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4858        @type where: C{None} or L{escript.FunctionSpace}
4859        @return: integral of arg.
4860        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4861      """      """
4862      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):  
4863         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4864      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4865         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
# Line 3654  def integrate(arg,where=None): Line 4870  def integrate(arg,where=None):
4870      else:      else:
4871        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4872    
4873    class Integrate_Symbol(DependendSymbol):
4874       """
4875       L{Symbol} representing the result of the spatial integration operator
4876       """
4877       def __init__(self,arg,where=None):
4878          """
4879          initialization of integration L{Symbol} with argument arg
4880          @param arg: argument of the integration
4881          @type arg: L{Symbol}.
4882          @param where: FunctionSpace in which the integration will be calculated.
4883                      If not present or C{None} an appropriate default is used.
4884          @type where: C{None} or L{escript.FunctionSpace}
4885          """
4886          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4887    
4888       def getMyCode(self,argstrs,format="escript"):
4889          """
4890          returns a program code that can be used to evaluate the symbol.
4891    
4892          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4893          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4894          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4895          @type format: C{str}
4896          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4897          @rtype: C{str}
4898          @raise NotImplementedError: if the requested format is not available
4899          """
4900          if format=="escript" or format=="str"  or format=="text":
4901             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4902          else:
4903             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4904    
4905       def substitute(self,argvals):
4906          """
4907          assigns new values to symbols in the definition of the symbol.
4908          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4909    
4910          @param argvals: new values assigned to symbols
4911          @type argvals: C{dict} with keywords of type L{Symbol}.
4912          @return: result of the substitution process. Operations are executed as much as possible.
4913          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4914          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4915          """
4916          if argvals.has_key(self):
4917             arg=argvals[self]
4918             if self.isAppropriateValue(arg):
4919                return arg
4920             else:
4921                raise TypeError,"%s: new value is not appropriate."%str(self)
4922          else:
4923             arg=self.getSubstitutedArguments(argvals)
4924             return integrate(arg[0],where=arg[1])
4925    
4926       def diff(self,arg):
4927          """
4928          differential of this object
4929    
4930          @param arg: the derivative is calculated with respect to arg
4931          @type arg: L{escript.Symbol}
4932          @return: derivative with respect to C{arg}
4933          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4934          """
4935          if arg==self:
4936             return identity(self.getShape())
4937          else:
4938             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4939    
4940    
4941  def interpolate(arg,where):  def interpolate(arg,where):
4942      """      """
4943      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
4944    
4945      @param arg:    interpolant      @param arg: interpolant
4946      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
4947        @param where: FunctionSpace to be interpolated to
4948        @type where: L{escript.FunctionSpace}
4949        @return: interpolated argument
4950        @rtype: C{escript.Data} or L{Symbol}
4951      """      """
4952      if testForZero(arg):      if isinstance(arg,Symbol):
4953        return 0         return Interpolate_Symbol(arg,where)
     elif isinstance(arg,Symbol):  
        return Interpolated_Symbol(arg,where)  
4954      else:      else:
4955         return escript.Data(arg,where)         return escript.Data(arg,where)
4956    
4957  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
4958      """     """
4959      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
4960       """
4961       def __init__(self,arg,where):
4962          """
4963          initialization of interpolation L{Symbol} with argument arg
4964          @param arg: argument of the interpolation
4965          @type arg: L{Symbol}.
4966          @param where: FunctionSpace into which the argument is interpolated.
4967          @type where: L{escript.FunctionSpace}
4968          """
4969          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4970    
4971      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
4972                    be calculated.        """
4973      @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)  
4974    
4975  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
4976      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4977      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.
4978          @type format: C{str}
4979          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4980          @rtype: C{str}
4981          @raise NotImplementedError: if the requested format is not available
4982          """
4983          if format=="escript" or format=="str"  or format=="text":
4984             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
4985          else:
4986             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4987    
4988      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
4989                    to be calculated.        """
4990      """        assigns new values to symbols in the definition of the symbol.
4991      d=arg.getDomain()        The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
     return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())  
4992    
4993  #=============================        @param argvals: new values assigned to symbols
4994  #        @type argvals: C{dict} with keywords of type L{Symbol}.
4995  # 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.
4996  # 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
4997  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4998          """
4999          if argvals.has_key(self):
5000             arg=argvals[self]
5001             if self.isAppropriateValue(arg):
5002                return arg
5003             else:
5004                raise TypeError,"%s: new value is not appropriate."%str(self)
5005          else:
5006             arg=self.getSubstitutedArguments(argvals)
5007             return interpolate(arg[0],where=arg[1])
5008    
5009       def diff(self,arg):
5010          """
5011          differential of this object
5012    
5013          @param arg: the derivative is calculated with respect to arg
5014          @type arg: L{escript.Symbol}
5015          @return: derivative with respect to C{arg}
5016          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
5017          """
5018          if arg==self:
5019             return identity(self.getShape())
5020          else:
5021             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5022    
 # functions involving the underlying Domain:  
5023    
5024  def transpose(arg,axis=None):  def div(arg,where=None):
5025      """      """
5026      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
5027    
5028      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
5029        @type arg: L{escript.Data} or L{Symbol}
5030        @param where: FunctionSpace in which the divergence will be calculated.
5031                      If not present or C{None} an appropriate default is used.
5032        @type where: C{None} or L{escript.FunctionSpace}
5033        @return: divergence of arg.
5034        @rtype: L{escript.Data} or L{Symbol}
5035      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
5036      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5037         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
5038      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
5039         # 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)  
5040      else:      else:
5041         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
5042        if not arg.getShape()==(dim,):
5043          raise ValueError,"div: expected shape is (%s,)"%dim
5044        return trace(grad(arg,where))
5045    
5046  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
5047      """      """
5048      Return      returns the jump of arg across the continuity of the domain
5049    
5050      @param arg:      @param arg: argument
5051        @type arg: L{escript.Data} or L{Symbol}
5052        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5053                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5054        @type domain: C{None} or L{escript.Domain}
5055        @return: jump of arg
5056        @rtype: L{escript.Data} or L{Symbol}
5057      """      """
5058      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
5059         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)  
5060    
5061    def L2(arg):
5062        """
5063        returns the L2 norm of arg at where
5064        
5065        @param arg: function which L2 to be calculated.
5066        @type arg: L{escript.Data} or L{Symbol}
5067        @return: L2 norm of arg.
5068        @rtype: L{float} or L{Symbol}
5069        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5070        """
5071        return sqrt(integrate(inner(arg,arg)))
5072    #=============================
5073    #
5074    
5075  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5076      """      """
5077      resorts the component of arg according to index      resorts the component of arg according to index
5078    
5079      """      """
5080      pass      raise NotImplementedError
5081  #  #
5082  # $Log: util.py,v $  # $Log: util.py,v $
5083  # Revision 1.14.2.16  2005/10/19 06:09:57  gross  # Revision 1.14.2.16  2005/10/19 06:09:57  gross

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

  ViewVC Help
Powered by ViewVC 1.1.26