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

  ViewVC Help
Powered by ViewVC 1.1.26