/[escript]/trunk-mpi-branch/escript/py_src/util.py
ViewVC logotype

Diff of /trunk-mpi-branch/escript/py_src/util.py

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

trunk/escript/py_src/util.py revision 341 by gross, Mon Dec 12 05:26:10 2005 UTC trunk-mpi-branch/escript/py_src/util.py revision 1295 by ksteube, Mon Sep 10 06:07:09 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: 329 $"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
23    
24    
# Line 32  import math Line 26  import math
26  import numarray  import numarray
27  import escript  import escript
28  import os  import os
29    from esys.escript import C_GeneralTensorProduct
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 maximum(arg0,arg1):  
 # def minimum(arg0,arg1):  
   
 # def transpose(arg,axis=None):  
 # def trace(arg,axis0=0,axis1=1):  
 # def reorderComponents(arg,index):  
   
 # def integrate(arg,where=None):  
 # def interpolate(arg,where):  
 # def div(arg,where=None):  
 # def grad(arg,where=None):  
   
 #  
 # slicing: get  
 #          set  
 #  
 # and derivatives  
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  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
75      """      """
76      writes a L{Data} objects into a files using the the VTK XML file format.      writes a L{Data} objects into a files using the the VTK XML file format.
77    
78      Example:      Example::
79    
80         tmp=Scalar(..)         tmp=Scalar(..)
81         v=Vector(..)         v=Vector(..)
# Line 84  def saveVTK(filename,domain=None,**data) Line 91  def saveVTK(filename,domain=None,**data)
91      @type <name>: L{Data} object.      @type <name>: L{Data} object.
92      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
93      """      """
94      if domain==None:      new_data={}
95         for i in data.keys():      for n,d in data.items():
96            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
97                fs=d.getFunctionSpace()
98                domain2=fs.getDomain()
99                if fs == escript.Solution(domain2):
100                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
101                elif fs == escript.ReducedSolution(domain2):
102                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
103                else:
104                   new_data[n]=d
105                if domain==None: domain=domain2
106      if domain==None:      if domain==None:
107          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
108      else:      domain.saveVTK(filename,new_data)
         domain.saveVTK(filename,data)  
109    
110  def saveDX(filename,domain=None,**data):  def saveDX(filename,domain=None,**data):
111      """      """
112      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.
113    
114      Example:      Example::
115    
116         tmp=Scalar(..)         tmp=Scalar(..)
117         v=Vector(..)         v=Vector(..)
# Line 112  def saveDX(filename,domain=None,**data): Line 127  def saveDX(filename,domain=None,**data):
127      @type <name>: L{Data} object.      @type <name>: L{Data} object.
128      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
129      """      """
130      if domain==None:      new_data={}
131         for i in data.keys():      for n,d in data.items():
132            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
133                fs=d.getFunctionSpace()
134                domain2=fs.getDomain()
135                if fs == escript.Solution(domain2):
136                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
137                elif fs == escript.ReducedSolution(domain2):
138                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
139                elif fs == escript.ContinuousFunction(domain2):
140                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
141                else:
142                   new_data[n]=d
143                if domain==None: domain=domain2
144      if domain==None:      if domain==None:
145          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
146      else:      domain.saveDX(filename,new_data)
         domain.saveDX(filename,data)  
147    
148  def kronecker(d=3):  def kronecker(d=3):
149     """     """
150     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
151    
152     @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
153     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
154     @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
155     @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}  
156     """     """
157     return identityTensor(d)     return identityTensor(d)
158    
# Line 143  def identity(shape=()): Line 167  def identity(shape=()):
167     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
168     """     """
169     if len(shape)>0:     if len(shape)>0:
170        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
171        if len(shape)==1:        if len(shape)==1:
172            for i0 in range(shape[0]):            for i0 in range(shape[0]):
173               out[i0,i0]=1.               out[i0,i0]=1.
   
174        elif len(shape)==2:        elif len(shape)==2:
175            for i0 in range(shape[0]):            for i0 in range(shape[0]):
176               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 163  def identityTensor(d=3): Line 186  def identityTensor(d=3):
186     return the dxd identity matrix     return the dxd identity matrix
187    
188     @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
189     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
190     @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
191     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
192     """     """
193     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
194        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
195     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
196           return identity((d.getDim(),))
197       else:
198           return identity((d,))
199    
200  def identityTensor4(d=3):  def identityTensor4(d=3):
201     """     """
# Line 178  def identityTensor4(d=3): Line 204  def identityTensor4(d=3):
204     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
205     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
206     @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
207     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
208     """     """
209     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
210        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
211     return identity((d,d))     elif isinstance(d,escript.Domain):
212           return identity((d.getDim(),d.getDim()))
213       else:
214           return identity((d,d))
215    
216  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
217     """     """
# Line 191  def unitVector(i=0,d=3): Line 220  def unitVector(i=0,d=3):
220     @param i: index     @param i: index
221     @type i: C{int}     @type i: C{int}
222     @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
223     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
224     @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
225     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
226     """     """
227     return kronecker(d)[i]     return kronecker(d)[i]
228    
# Line 249  def inf(arg): Line 278  def inf(arg):
278    
279      @param arg: argument      @param arg: argument
280      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
281      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
282      @rtype: C{float}      @rtype: C{float}
283      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
284      """      """
# Line 268  def inf(arg): Line 297  def inf(arg):
297  #=========================================================================  #=========================================================================
298  #   some little helpers  #   some little helpers
299  #=========================================================================  #=========================================================================
300  def pokeShape(arg):  def getRank(arg):
301        """
302        identifies the rank of its argument
303    
304        @param arg: a given object
305        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
306        @return: the rank of the argument
307        @rtype: C{int}
308        @raise TypeError: if type of arg cannot be processed
309        """
310    
311        if isinstance(arg,numarray.NumArray):
312            return arg.rank
313        elif isinstance(arg,escript.Data):
314            return arg.getRank()
315        elif isinstance(arg,float):
316            return 0
317        elif isinstance(arg,int):
318            return 0
319        elif isinstance(arg,Symbol):
320            return arg.getRank()
321        else:
322          raise TypeError,"getShape: cannot identify shape"
323    def getShape(arg):
324      """      """
325      identifies the shape of its argument      identifies the shape of its argument
326    
# Line 290  def pokeShape(arg): Line 342  def pokeShape(arg):
342      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
343          return arg.getShape()          return arg.getShape()
344      else:      else:
345        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
346    
347  def pokeDim(arg):  def pokeDim(arg):
348      """      """
# Line 313  def commonShape(arg0,arg1): Line 365  def commonShape(arg0,arg1):
365      """      """
366      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.
367    
368      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
369      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
370      @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.
371      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
372      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
373      """      """
374      sh0=pokeShape(arg0)      sh0=getShape(arg0)
375      sh1=pokeShape(arg1)      sh1=getShape(arg1)
376      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
377         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
378               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
# Line 338  def commonDim(*args): Line 390  def commonDim(*args):
390      """      """
391      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.
392    
393      @param *args: given objects      @param args: given objects
394      @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
395               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
396      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 360  def testForZero(arg): Line 412  def testForZero(arg):
412    
413      @param arg: a given object      @param arg: a given object
414      @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}
415      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
416      @rtype : C{bool}      @rtype: C{bool}
417      """      """
418      try:      if isinstance(arg,numarray.NumArray):
419           return not Lsup(arg)>0.
420        elif isinstance(arg,escript.Data):
421           return False
422        elif isinstance(arg,float):
423           return not Lsup(arg)>0.
424        elif isinstance(arg,int):
425         return not Lsup(arg)>0.         return not Lsup(arg)>0.
426      except TypeError:      elif isinstance(arg,Symbol):
427           return False
428        else:
429         return False         return False
430    
431  def matchType(arg0=0.,arg1=0.):  def matchType(arg0=0.,arg1=0.):
# Line 386  def matchType(arg0=0.,arg1=0.): Line 446  def matchType(arg0=0.,arg1=0.):
446         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
447            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
448         elif isinstance(arg1,float):         elif isinstance(arg1,float):
449            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
450         elif isinstance(arg1,int):         elif isinstance(arg1,int):
451            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
452         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
453            pass            pass
454         else:         else:
# Line 412  def matchType(arg0=0.,arg1=0.): Line 472  def matchType(arg0=0.,arg1=0.):
472         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
473            pass            pass
474         elif isinstance(arg1,float):         elif isinstance(arg1,float):
475            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
476         elif isinstance(arg1,int):         elif isinstance(arg1,int):
477            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
478         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
479            pass            pass
480         else:         else:
481            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
482      elif isinstance(arg0,float):      elif isinstance(arg0,float):
483         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
484            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
485         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
486            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
487         elif isinstance(arg1,float):         elif isinstance(arg1,float):
488            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
489            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
490         elif isinstance(arg1,int):         elif isinstance(arg1,int):
491            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
492            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
493         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
494            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
495         else:         else:
496            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
497      elif isinstance(arg0,int):      elif isinstance(arg0,int):
498         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
499            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
500         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
501            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
502         elif isinstance(arg1,float):         elif isinstance(arg1,float):
503            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
504            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
505         elif isinstance(arg1,int):         elif isinstance(arg1,int):
506            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
507            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
508         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
509            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
510         else:         else:
511            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
512      else:      else:
# Line 456  def matchType(arg0=0.,arg1=0.): Line 516  def matchType(arg0=0.,arg1=0.):
516    
517  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
518      """      """
519            return representations of arg0 amd arg1 which ahve the same shape
520    
521      If shape is not given the shape "largest" shape of args is used.      @param arg0: a given object
522        @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
523      @param args: a given ob      @param arg1: a given object
524      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
525      @return: True if the argument is identical to zero.      @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
526      @rtype: C{list} of C{int}      @rtype: C{tuple}
527      """      """
528      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
529      sh0=pokeShape(arg0)      sh0=getShape(arg0)
530      sh1=pokeShape(arg1)      sh1=getShape(arg1)
531      if len(sh0)<len(sh):      if len(sh0)<len(sh):
532         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
533      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
534         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
535      else:      else:
536         return arg0,arg1         return arg0,arg1
537  #=========================================================  #=========================================================
# Line 491  class Symbol(object): Line 551  class Symbol(object):
551         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
552         symbols or any other object.         symbols or any other object.
553    
554         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
555         @type arg: C{list}         @type args: C{list}
556         @param shape: the shape         @param shape: the shape
557         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
558         @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
# Line 535  class Symbol(object): Line 595  class Symbol(object):
595         """         """
596         the shape of the symbol.         the shape of the symbol.
597    
598         @return : the shape of the symbol.         @return: the shape of the symbol.
599         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
600         """         """
601         return self.__shape         return self.__shape
# Line 544  class Symbol(object): Line 604  class Symbol(object):
604         """         """
605         the spatial dimension         the spatial dimension
606    
607         @return : the spatial dimension         @return: the spatial dimension
608         @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.
609         """         """
610         return self.__dim         return self.__dim
# Line 568  class Symbol(object): Line 628  class Symbol(object):
628         """         """
629         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.
630    
631         @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].
632         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
633         @rtype: C{list} of objects         @rtype: C{list} of objects
634         @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.
635         """         """
636         out=[]         out=[]
637         for a in self.getArgument():         for a in self.getArgument():
# Line 595  class Symbol(object): Line 655  class Symbol(object):
655            if isinstance(a,Symbol):            if isinstance(a,Symbol):
656               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
657            else:            else:
658                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
659                if len(s)>0:                if len(s)>0:
660                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
661                else:                else:
662                   out.append(a)                   out.append(a)
663         return out         return out
# Line 687  class Symbol(object): Line 747  class Symbol(object):
747         else:         else:
748            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
749            if len(s)>0:            if len(s)>0:
750               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
751            else:            else:
752               return 0.               return 0.
753    
# Line 695  class Symbol(object): Line 755  class Symbol(object):
755         """         """
756         returns -self.         returns -self.
757    
758         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
759         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
760         """         """
761         return self*(-1.)         return self*(-1.)
# Line 704  class Symbol(object): Line 764  class Symbol(object):
764         """         """
765         returns +self.         returns +self.
766    
767         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
768         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
769         """         """
770         return self*(1.)         return self*(1.)
771    
772     def __abs__(self):     def __abs__(self):
773         """         """
774         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
775         """         """
776         return Abs_Symbol(self)         return Abs_Symbol(self)
777    
# Line 721  class Symbol(object): Line 781  class Symbol(object):
781    
782         @param other: object to be added to this object         @param other: object to be added to this object
783         @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}.
784         @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}
785         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
786         """         """
787         return add(self,other)         return add(self,other)
# Line 732  class Symbol(object): Line 792  class Symbol(object):
792    
793         @param other: object this object is added to         @param other: object this object is added to
794         @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}.
795         @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
796         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
797         """         """
798         return add(other,self)         return add(other,self)
# Line 743  class Symbol(object): Line 803  class Symbol(object):
803    
804         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
805         @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}.
806         @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
807         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
808         """         """
809         return add(self,-other)         return add(self,-other)
# Line 754  class Symbol(object): Line 814  class Symbol(object):
814    
815         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
816         @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}.
817         @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}.
818         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
819         """         """
820         return add(-self,other)         return add(-self,other)
# Line 765  class Symbol(object): Line 825  class Symbol(object):
825    
826         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
827         @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}.
828         @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}.
829         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
830         """         """
831         return mult(self,other)         return mult(self,other)
# Line 776  class Symbol(object): Line 836  class Symbol(object):
836    
837         @param other: object this object is multiplied with         @param other: object this object is multiplied with
838         @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}.
839         @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.
840         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
841         """         """
842         return mult(other,self)         return mult(other,self)
# Line 787  class Symbol(object): Line 847  class Symbol(object):
847    
848         @param other: object dividing this object         @param other: object dividing this object
849         @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}.
850         @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}
851         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
852         """         """
853         return quotient(self,other)         return quotient(self,other)
# Line 798  class Symbol(object): Line 858  class Symbol(object):
858    
859         @param other: object dividing this object         @param other: object dividing this object
860         @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}.
861         @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
862         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
863         """         """
864         return quotient(other,self)         return quotient(other,self)
# Line 809  class Symbol(object): Line 869  class Symbol(object):
869    
870         @param other: exponent         @param other: exponent
871         @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}.
872         @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}
873         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
874         """         """
875         return power(self,other)         return power(self,other)
# Line 820  class Symbol(object): Line 880  class Symbol(object):
880    
881         @param other: basis         @param other: basis
882         @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}.
883         @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
884         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
885         """         """
886         return power(other,self)         return power(other,self)
887    
888       def __getitem__(self,index):
889           """
890           returns the slice defined by index
891    
892           @param index: defines a
893           @type index: C{slice} or C{int} or a C{tuple} of them
894           @return: a L{Symbol} representing the slice defined by index
895           @rtype: L{DependendSymbol}
896           """
897           return GetSlice_Symbol(self,index)
898    
899  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
900     """     """
901     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.
902     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  
903        
904     Example:     Example::
905        
906     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
907     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
908     print u1==u2       print u1==u2
909     False       False
910        
911        but       but::
912    
913     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
914     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
915     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
916     print u1==u2, u1==u3       print u1==u2, u1==u3
917     True False       True False
918    
919     @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.
920     """     """
# Line 875  class DependendSymbol(Symbol): Line 946  class DependendSymbol(Symbol):
946  #=========================================================  #=========================================================
947  #  Unary operations prserving the shape  #  Unary operations prserving the shape
948  #========================================================  #========================================================
949    class GetSlice_Symbol(DependendSymbol):
950       """
951       L{Symbol} representing getting a slice for a L{Symbol}
952       """
953       def __init__(self,arg,index):
954          """
955          initialization of wherePositive L{Symbol} with argument arg
956          @param arg: argument
957          @type arg: L{Symbol}.
958          @param index: defines index
959          @type index: C{slice} or C{int} or a C{tuple} of them
960          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
961          @raises ValueError: if a step is given
962          """
963          if not isinstance(index,tuple): index=(index,)
964          if len(index)>arg.getRank():
965               raise IndexError,"GetSlice_Symbol: index out of range."
966          sh=()
967          index2=()
968          for i in range(len(index)):
969             ix=index[i]
970             if isinstance(ix,int):
971                if ix<0 or ix>=arg.getShape()[i]:
972                   raise ValueError,"GetSlice_Symbol: index out of range."
973                index2=index2+(ix,)
974             else:
975               if not ix.step==None:
976                 raise ValueError,"GetSlice_Symbol: steping is not supported."
977               if ix.start==None:
978                  s=0
979               else:
980                  s=ix.start
981               if ix.stop==None:
982                  e=arg.getShape()[i]
983               else:
984                  e=ix.stop
985                  if e>arg.getShape()[i]:
986                     raise IndexError,"GetSlice_Symbol: index out of range."
987               index2=index2+(slice(s,e),)
988               if e>s:
989                   sh=sh+(e-s,)
990               elif s>e:
991                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
992          for i in range(len(index),arg.getRank()):
993              index2=index2+(slice(0,arg.getShape()[i]),)
994              sh=sh+(arg.getShape()[i],)
995          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
996    
997       def getMyCode(self,argstrs,format="escript"):
998          """
999          returns a program code that can be used to evaluate the symbol.
1000    
1001          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1002          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1003          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1004          @type format: C{str}
1005          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1006          @rtype: C{str}
1007          @raise NotImplementedError: if the requested format is not available
1008          """
1009          if format=="escript" or format=="str"  or format=="text":
1010             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
1011          else:
1012             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
1013    
1014       def substitute(self,argvals):
1015          """
1016          assigns new values to symbols in the definition of the symbol.
1017          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1018    
1019          @param argvals: new values assigned to symbols
1020          @type argvals: C{dict} with keywords of type L{Symbol}.
1021          @return: result of the substitution process. Operations are executed as much as possible.
1022          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1023          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1024          """
1025          if argvals.has_key(self):
1026             arg=argvals[self]
1027             if self.isAppropriateValue(arg):
1028                return arg
1029             else:
1030                raise TypeError,"%s: new value is not appropriate."%str(self)
1031          else:
1032             args=self.getSubstitutedArguments(argvals)
1033             arg=args[0]
1034             index=args[1]
1035             return arg.__getitem__(index)
1036    
1037  def log10(arg):  def log10(arg):
1038     """     """
1039     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
1040    
1041     @param arg: argument     @param arg: argument
1042     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1043     @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.
1044     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1045     """     """
1046     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 903  def wherePositive(arg): Line 1062  def wherePositive(arg):
1062    
1063     @param arg: argument     @param arg: argument
1064     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1065     @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.
1066     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1067     """     """
1068     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1069        if arg.rank==0:        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1070           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1071             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))  
1072     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1073        return arg._wherePositive()        return arg._wherePositive()
1074     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 953  class WherePositive_Symbol(DependendSymb Line 1108  class WherePositive_Symbol(DependendSymb
1108        @type format: C{str}        @type format: C{str}
1109        @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.
1110        @rtype: C{str}        @rtype: C{str}
1111        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1112        """        """
1113        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1114            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 989  def whereNegative(arg): Line 1144  def whereNegative(arg):
1144    
1145     @param arg: argument     @param arg: argument
1146     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1147     @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.
1148     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1149     """     """
1150     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1151        if arg.rank==0:        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1152           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1153             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))  
1154     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1155        return arg._whereNegative()        return arg._whereNegative()
1156     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1039  class WhereNegative_Symbol(DependendSymb Line 1190  class WhereNegative_Symbol(DependendSymb
1190        @type format: C{str}        @type format: C{str}
1191        @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.
1192        @rtype: C{str}        @rtype: C{str}
1193        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1194        """        """
1195        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1196            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1075  def whereNonNegative(arg): Line 1226  def whereNonNegative(arg):
1226    
1227     @param arg: argument     @param arg: argument
1228     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1229     @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.
1230     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1231     """     """
1232     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1233        if arg.rank==0:        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1234           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1235             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))  
1236     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1237        return arg._whereNonNegative()        return arg._whereNonNegative()
1238     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1109  def whereNonPositive(arg): Line 1256  def whereNonPositive(arg):
1256    
1257     @param arg: argument     @param arg: argument
1258     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1259     @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.
1260     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1261     """     """
1262     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1263        if arg.rank==0:        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1264           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1265             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.  
1266     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1267        return arg._whereNonPositive()        return arg._whereNonPositive()
1268     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1145  def whereZero(arg,tol=0.): Line 1288  def whereZero(arg,tol=0.):
1288     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1289     @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.
1290     @type tol: C{float}     @type tol: C{float}
1291     @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.
1292     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1293     """     """
1294     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1295        if arg.rank==0:        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1296           if abs(arg)<=tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1297             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
1298     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1299        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1300     elif isinstance(arg,float):     elif isinstance(arg,float):
1301        if abs(arg)<=tol:        if abs(arg)<=tol:
1302          return 1.          return 1.
# Line 1198  class WhereZero_Symbol(DependendSymbol): Line 1334  class WhereZero_Symbol(DependendSymbol):
1334        @type format: C{str}        @type format: C{str}
1335        @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.
1336        @rtype: C{str}        @rtype: C{str}
1337        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1338        """        """
1339        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1340           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1232  def whereNonZero(arg,tol=0.): Line 1368  def whereNonZero(arg,tol=0.):
1368    
1369     @param arg: argument     @param arg: argument
1370     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1371     @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.
1372     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1373     """     """
1374     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1375        if arg.rank==0:        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1376          if abs(arg)>tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1377             return numarray.array(1.)        return out
         else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
1378     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1379        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1380     elif isinstance(arg,float):     elif isinstance(arg,float):
1381        if abs(arg)>tol:        if abs(arg)>tol:
1382          return 1.          return 1.
# Line 1263  def whereNonZero(arg,tol=0.): Line 1392  def whereNonZero(arg,tol=0.):
1392     else:     else:
1393        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1394    
1395    def erf(arg):
1396       """
1397       returns erf of argument arg
1398    
1399       @param arg: argument
1400       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1401       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1402       @raises TypeError: if the type of the argument is not expected.
1403       """
1404       if isinstance(arg,escript.Data):
1405          return arg._erf()
1406       else:
1407          raise TypeError,"erf: Unknown argument type."
1408    
1409  def sin(arg):  def sin(arg):
1410     """     """
1411     returns sine of argument arg     returns sine of argument arg
1412    
1413     @param arg: argument     @param arg: argument
1414     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1415     @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.
1416     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1417     """     """
1418     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1307  class Sin_Symbol(DependendSymbol): Line 1450  class Sin_Symbol(DependendSymbol):
1450        @type format: C{str}        @type format: C{str}
1451        @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.
1452        @rtype: C{str}        @rtype: C{str}
1453        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1454        """        """
1455        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1456            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1359  def cos(arg): Line 1502  def cos(arg):
1502    
1503     @param arg: argument     @param arg: argument
1504     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1505     @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.
1506     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1507     """     """
1508     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1397  class Cos_Symbol(DependendSymbol): Line 1540  class Cos_Symbol(DependendSymbol):
1540        @type format: C{str}        @type format: C{str}
1541        @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.
1542        @rtype: C{str}        @rtype: C{str}
1543        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1544        """        """
1545        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1546            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1449  def tan(arg): Line 1592  def tan(arg):
1592    
1593     @param arg: argument     @param arg: argument
1594     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1595     @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.
1596     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1597     """     """
1598     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1487  class Tan_Symbol(DependendSymbol): Line 1630  class Tan_Symbol(DependendSymbol):
1630        @type format: C{str}        @type format: C{str}
1631        @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.
1632        @rtype: C{str}        @rtype: C{str}
1633        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1634        """        """
1635        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1636            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1539  def asin(arg): Line 1682  def asin(arg):
1682    
1683     @param arg: argument     @param arg: argument
1684     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1685     @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.
1686     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1687     """     """
1688     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1577  class Asin_Symbol(DependendSymbol): Line 1720  class Asin_Symbol(DependendSymbol):
1720        @type format: C{str}        @type format: C{str}
1721        @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.
1722        @rtype: C{str}        @rtype: C{str}
1723        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1724        """        """
1725        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1726            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1629  def acos(arg): Line 1772  def acos(arg):
1772    
1773     @param arg: argument     @param arg: argument
1774     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1775     @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.
1776     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1777     """     """
1778     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1667  class Acos_Symbol(DependendSymbol): Line 1810  class Acos_Symbol(DependendSymbol):
1810        @type format: C{str}        @type format: C{str}
1811        @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.
1812        @rtype: C{str}        @rtype: C{str}
1813        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1814        """        """
1815        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1816            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1719  def atan(arg): Line 1862  def atan(arg):
1862    
1863     @param arg: argument     @param arg: argument
1864     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1865     @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.
1866     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1867     """     """
1868     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1757  class Atan_Symbol(DependendSymbol): Line 1900  class Atan_Symbol(DependendSymbol):
1900        @type format: C{str}        @type format: C{str}
1901        @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.
1902        @rtype: C{str}        @rtype: C{str}
1903        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1904        """        """
1905        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1906            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1809  def sinh(arg): Line 1952  def sinh(arg):
1952    
1953     @param arg: argument     @param arg: argument
1954     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1955     @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.
1956     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1957     """     """
1958     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1847  class Sinh_Symbol(DependendSymbol): Line 1990  class Sinh_Symbol(DependendSymbol):
1990        @type format: C{str}        @type format: C{str}
1991        @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.
1992        @rtype: C{str}        @rtype: C{str}
1993        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1994        """        """
1995        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1996            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1899  def cosh(arg): Line 2042  def cosh(arg):
2042    
2043     @param arg: argument     @param arg: argument
2044     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2045     @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.
2046     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2047     """     """
2048     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1937  class Cosh_Symbol(DependendSymbol): Line 2080  class Cosh_Symbol(DependendSymbol):
2080        @type format: C{str}        @type format: C{str}
2081        @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.
2082        @rtype: C{str}        @rtype: C{str}
2083        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2084        """        """
2085        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2086            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1989  def tanh(arg): Line 2132  def tanh(arg):
2132    
2133     @param arg: argument     @param arg: argument
2134     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2135     @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.
2136     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2137     """     """
2138     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2027  class Tanh_Symbol(DependendSymbol): Line 2170  class Tanh_Symbol(DependendSymbol):
2170        @type format: C{str}        @type format: C{str}
2171        @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.
2172        @rtype: C{str}        @rtype: C{str}
2173        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2174        """        """
2175        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2176            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2079  def asinh(arg): Line 2222  def asinh(arg):
2222    
2223     @param arg: argument     @param arg: argument
2224     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2225     @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.
2226     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2227     """     """
2228     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2117  class Asinh_Symbol(DependendSymbol): Line 2260  class Asinh_Symbol(DependendSymbol):
2260        @type format: C{str}        @type format: C{str}
2261        @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.
2262        @rtype: C{str}        @rtype: C{str}
2263        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2264        """        """
2265        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2266            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2169  def acosh(arg): Line 2312  def acosh(arg):
2312    
2313     @param arg: argument     @param arg: argument
2314     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2315     @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.
2316     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2317     """     """
2318     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2207  class Acosh_Symbol(DependendSymbol): Line 2350  class Acosh_Symbol(DependendSymbol):
2350        @type format: C{str}        @type format: C{str}
2351        @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.
2352        @rtype: C{str}        @rtype: C{str}
2353        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2354        """        """
2355        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2356            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2259  def atanh(arg): Line 2402  def atanh(arg):
2402    
2403     @param arg: argument     @param arg: argument
2404     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2405     @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.
2406     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2407     """     """
2408     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2297  class Atanh_Symbol(DependendSymbol): Line 2440  class Atanh_Symbol(DependendSymbol):
2440        @type format: C{str}        @type format: C{str}
2441        @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.
2442        @rtype: C{str}        @rtype: C{str}
2443        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2444        """        """
2445        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2446            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2349  def exp(arg): Line 2492  def exp(arg):
2492    
2493     @param arg: argument     @param arg: argument
2494     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2495     @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.
2496     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2497     """     """
2498     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2387  class Exp_Symbol(DependendSymbol): Line 2530  class Exp_Symbol(DependendSymbol):
2530        @type format: C{str}        @type format: C{str}
2531        @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.
2532        @rtype: C{str}        @rtype: C{str}
2533        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2534        """        """
2535        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2536            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2439  def sqrt(arg): Line 2582  def sqrt(arg):
2582    
2583     @param arg: argument     @param arg: argument
2584     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2585     @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.
2586     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2587     """     """
2588     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2477  class Sqrt_Symbol(DependendSymbol): Line 2620  class Sqrt_Symbol(DependendSymbol):
2620        @type format: C{str}        @type format: C{str}
2621        @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.
2622        @rtype: C{str}        @rtype: C{str}
2623        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2624        """        """
2625        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2626            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2529  def log(arg): Line 2672  def log(arg):
2672    
2673     @param arg: argument     @param arg: argument
2674     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2675     @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.
2676     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2677     """     """
2678     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2567  class Log_Symbol(DependendSymbol): Line 2710  class Log_Symbol(DependendSymbol):
2710        @type format: C{str}        @type format: C{str}
2711        @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.
2712        @rtype: C{str}        @rtype: C{str}
2713        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2714        """        """
2715        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2716            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2619  def sign(arg): Line 2762  def sign(arg):
2762    
2763     @param arg: argument     @param arg: argument
2764     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2765     @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.
2766     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2767     """     """
2768     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2667  class Abs_Symbol(DependendSymbol): Line 2810  class Abs_Symbol(DependendSymbol):
2810        @type format: C{str}        @type format: C{str}
2811        @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.
2812        @rtype: C{str}        @rtype: C{str}
2813        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2814        """        """
2815        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2816            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2719  def minval(arg): Line 2862  def minval(arg):
2862    
2863     @param arg: argument     @param arg: argument
2864     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2865     @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.
2866     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2867     """     """
2868     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2760  class Minval_Symbol(DependendSymbol): Line 2903  class Minval_Symbol(DependendSymbol):
2903        @type format: C{str}        @type format: C{str}
2904        @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.
2905        @rtype: C{str}        @rtype: C{str}
2906        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2907        """        """
2908        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2909            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2796  def maxval(arg): Line 2939  def maxval(arg):
2939    
2940     @param arg: argument     @param arg: argument
2941     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2942     @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.
2943     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2944     """     """
2945     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2837  class Maxval_Symbol(DependendSymbol): Line 2980  class Maxval_Symbol(DependendSymbol):
2980        @type format: C{str}        @type format: C{str}
2981        @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.
2982        @rtype: C{str}        @rtype: C{str}
2983        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2984        """        """
2985        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2986            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2873  def length(arg): Line 3016  def length(arg):
3016    
3017     @param arg: argument     @param arg: argument
3018     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3019     @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.
3020     """     """
3021     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
3022    
3023    def trace(arg,axis_offset=0):
3024       """
3025       returns the trace of arg which the sum of arg[k,k] over k.
3026    
3027       @param arg: argument
3028       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3029       @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
3030                      C{axis_offset} and axis_offset+1 must be equal.
3031       @type axis_offset: C{int}
3032       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
3033       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3034       """
3035       if isinstance(arg,numarray.NumArray):
3036          sh=arg.shape
3037          if len(sh)<2:
3038            raise ValueError,"rank of argument must be greater than 1"
3039          if axis_offset<0 or axis_offset>len(sh)-2:
3040            raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
3041          s1=1
3042          for i in range(axis_offset): s1*=sh[i]
3043          s2=1
3044          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3045          if not sh[axis_offset] == sh[axis_offset+1]:
3046            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3047          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
3048          out=numarray.zeros([s1,s2],numarray.Float64)
3049          for i1 in range(s1):
3050            for i2 in range(s2):
3051                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
3052          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3053          return out
3054       elif isinstance(arg,escript.Data):
3055          if arg.getRank()<2:
3056            raise ValueError,"rank of argument must be greater than 1"
3057          if axis_offset<0 or axis_offset>arg.getRank()-2:
3058            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3059          s=list(arg.getShape())        
3060          if not s[axis_offset] == s[axis_offset+1]:
3061            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3062          return arg._trace(axis_offset)
3063       elif isinstance(arg,float):
3064          raise TypeError,"illegal argument type float."
3065       elif isinstance(arg,int):
3066          raise TypeError,"illegal argument type int."
3067       elif isinstance(arg,Symbol):
3068          return Trace_Symbol(arg,axis_offset)
3069       else:
3070          raise TypeError,"Unknown argument type."
3071    
3072    class Trace_Symbol(DependendSymbol):
3073       """
3074       L{Symbol} representing the result of the trace function
3075       """
3076       def __init__(self,arg,axis_offset=0):
3077          """
3078          initialization of trace L{Symbol} with argument arg
3079          @param arg: argument of function
3080          @type arg: L{Symbol}.
3081          @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
3082                      C{axis_offset} and axis_offset+1 must be equal.
3083          @type axis_offset: C{int}
3084          """
3085          if arg.getRank()<2:
3086            raise ValueError,"rank of argument must be greater than 1"
3087          if axis_offset<0 or axis_offset>arg.getRank()-2:
3088            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3089          s=list(arg.getShape())        
3090          if not s[axis_offset] == s[axis_offset+1]:
3091            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3092          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3093    
3094       def getMyCode(self,argstrs,format="escript"):
3095          """
3096          returns a program code that can be used to evaluate the symbol.
3097    
3098          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3099          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3100          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3101          @type format: C{str}
3102          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3103          @rtype: C{str}
3104          @raise NotImplementedError: if the requested format is not available
3105          """
3106          if format=="escript" or format=="str"  or format=="text":
3107             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3108          else:
3109             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3110    
3111       def substitute(self,argvals):
3112          """
3113          assigns new values to symbols in the definition of the symbol.
3114          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3115    
3116          @param argvals: new values assigned to symbols
3117          @type argvals: C{dict} with keywords of type L{Symbol}.
3118          @return: result of the substitution process. Operations are executed as much as possible.
3119          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3120          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3121          """
3122          if argvals.has_key(self):
3123             arg=argvals[self]
3124             if self.isAppropriateValue(arg):
3125                return arg
3126             else:
3127                raise TypeError,"%s: new value is not appropriate."%str(self)
3128          else:
3129             arg=self.getSubstitutedArguments(argvals)
3130             return trace(arg[0],axis_offset=arg[1])
3131    
3132       def diff(self,arg):
3133          """
3134          differential of this object
3135    
3136          @param arg: the derivative is calculated with respect to arg
3137          @type arg: L{escript.Symbol}
3138          @return: derivative with respect to C{arg}
3139          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3140          """
3141          if arg==self:
3142             return identity(self.getShape())
3143          else:
3144             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3145    
3146    def transpose(arg,axis_offset=None):
3147       """
3148       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3149    
3150       @param arg: argument
3151       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3152       @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.
3153                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3154       @type axis_offset: C{int}
3155       @return: transpose of arg
3156       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3157       """
3158       if isinstance(arg,numarray.NumArray):
3159          if axis_offset==None: axis_offset=int(arg.rank/2)
3160          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3161       elif isinstance(arg,escript.Data):
3162          r=arg.getRank()
3163          if axis_offset==None: axis_offset=int(r/2)
3164          if axis_offset<0 or axis_offset>r:
3165            raise ValueError,"axis_offset must be between 0 and %s"%r
3166          return arg._transpose(axis_offset)
3167       elif isinstance(arg,float):
3168          if not ( axis_offset==0 or axis_offset==None):
3169            raise ValueError,"axis_offset must be 0 for float argument"
3170          return arg
3171       elif isinstance(arg,int):
3172          if not ( axis_offset==0 or axis_offset==None):
3173            raise ValueError,"axis_offset must be 0 for int argument"
3174          return float(arg)
3175       elif isinstance(arg,Symbol):
3176          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3177          return Transpose_Symbol(arg,axis_offset)
3178       else:
3179          raise TypeError,"Unknown argument type."
3180    
3181    class Transpose_Symbol(DependendSymbol):
3182       """
3183       L{Symbol} representing the result of the transpose function
3184       """
3185       def __init__(self,arg,axis_offset=None):
3186          """
3187          initialization of transpose L{Symbol} with argument arg
3188    
3189          @param arg: argument of function
3190          @type arg: L{Symbol}.
3191          @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.
3192                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3193          @type axis_offset: C{int}
3194          """
3195          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3196          if axis_offset<0 or axis_offset>arg.getRank():
3197            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3198          s=arg.getShape()
3199          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3200    
3201       def getMyCode(self,argstrs,format="escript"):
3202          """
3203          returns a program code that can be used to evaluate the symbol.
3204    
3205          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3206          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3207          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3208          @type format: C{str}
3209          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3210          @rtype: C{str}
3211          @raise NotImplementedError: if the requested format is not available
3212          """
3213          if format=="escript" or format=="str"  or format=="text":
3214             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3215          else:
3216             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3217    
3218       def substitute(self,argvals):
3219          """
3220          assigns new values to symbols in the definition of the symbol.
3221          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3222    
3223          @param argvals: new values assigned to symbols
3224          @type argvals: C{dict} with keywords of type L{Symbol}.
3225          @return: result of the substitution process. Operations are executed as much as possible.
3226          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3227          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3228          """
3229          if argvals.has_key(self):
3230             arg=argvals[self]
3231             if self.isAppropriateValue(arg):
3232                return arg
3233             else:
3234                raise TypeError,"%s: new value is not appropriate."%str(self)
3235          else:
3236             arg=self.getSubstitutedArguments(argvals)
3237             return transpose(arg[0],axis_offset=arg[1])
3238    
3239       def diff(self,arg):
3240          """
3241          differential of this object
3242    
3243          @param arg: the derivative is calculated with respect to arg
3244          @type arg: L{escript.Symbol}
3245          @return: derivative with respect to C{arg}
3246          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3247          """
3248          if arg==self:
3249             return identity(self.getShape())
3250          else:
3251             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3252    
3253    def swap_axes(arg,axis0=0,axis1=1):
3254       """
3255       returns the swap of arg by swaping the components axis0 and axis1
3256    
3257       @param arg: argument
3258       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3259       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3260       @type axis0: C{int}
3261       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3262       @type axis1: C{int}
3263       @return: C{arg} with swaped components
3264       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3265       """
3266       if axis0 > axis1:
3267          axis0,axis1=axis1,axis0
3268       if isinstance(arg,numarray.NumArray):
3269          return numarray.swapaxes(arg,axis0,axis1)
3270       elif isinstance(arg,escript.Data):
3271          return arg._swap_axes(axis0,axis1)
3272       elif isinstance(arg,float):
3273          raise TyepError,"float argument is not supported."
3274       elif isinstance(arg,int):
3275          raise TyepError,"int argument is not supported."
3276       elif isinstance(arg,Symbol):
3277          return SwapAxes_Symbol(arg,axis0,axis1)
3278       else:
3279          raise TypeError,"Unknown argument type."
3280    
3281    class SwapAxes_Symbol(DependendSymbol):
3282       """
3283       L{Symbol} representing the result of the swap function
3284       """
3285       def __init__(self,arg,axis0=0,axis1=1):
3286          """
3287          initialization of swap L{Symbol} with argument arg
3288    
3289          @param arg: argument
3290          @type arg: L{Symbol}.
3291          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3292          @type axis0: C{int}
3293          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3294          @type axis1: C{int}
3295          """
3296          if arg.getRank()<2:
3297             raise ValueError,"argument must have at least rank 2."
3298          if axis0<0 or axis0>arg.getRank()-1:
3299             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3300          if axis1<0 or axis1>arg.getRank()-1:
3301             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3302          if axis0 == axis1:
3303             raise ValueError,"axis indices must be different."
3304          if axis0 > axis1:
3305             axis0,axis1=axis1,axis0
3306          s=arg.getShape()
3307          s_out=[]
3308          for i in range(len(s)):
3309             if i == axis0:
3310                s_out.append(s[axis1])
3311             elif i == axis1:
3312                s_out.append(s[axis0])
3313             else:
3314                s_out.append(s[i])
3315          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3316    
3317       def getMyCode(self,argstrs,format="escript"):
3318          """
3319          returns a program code that can be used to evaluate the symbol.
3320    
3321          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3322          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3323          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3324          @type format: C{str}
3325          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3326          @rtype: C{str}
3327          @raise NotImplementedError: if the requested format is not available
3328          """
3329          if format=="escript" or format=="str"  or format=="text":
3330             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3331          else:
3332             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3333    
3334       def substitute(self,argvals):
3335          """
3336          assigns new values to symbols in the definition of the symbol.
3337          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3338    
3339          @param argvals: new values assigned to symbols
3340          @type argvals: C{dict} with keywords of type L{Symbol}.
3341          @return: result of the substitution process. Operations are executed as much as possible.
3342          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3343          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3344          """
3345          if argvals.has_key(self):
3346             arg=argvals[self]
3347             if self.isAppropriateValue(arg):
3348                return arg
3349             else:
3350                raise TypeError,"%s: new value is not appropriate."%str(self)
3351          else:
3352             arg=self.getSubstitutedArguments(argvals)
3353             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3354    
3355       def diff(self,arg):
3356          """
3357          differential of this object
3358    
3359          @param arg: the derivative is calculated with respect to arg
3360          @type arg: L{escript.Symbol}
3361          @return: derivative with respect to C{arg}
3362          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3363          """
3364          if arg==self:
3365             return identity(self.getShape())
3366          else:
3367             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3368    
3369    def symmetric(arg):
3370        """
3371        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3372    
3373        @param arg: square matrix. Must have rank 2 or 4 and be square.
3374        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3375        @return: symmetric part of arg
3376        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3377        """
3378        if isinstance(arg,numarray.NumArray):
3379          if arg.rank==2:
3380            if not (arg.shape[0]==arg.shape[1]):
3381               raise ValueError,"argument must be square."
3382          elif arg.rank==4:
3383            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3384               raise ValueError,"argument must be square."
3385          else:
3386            raise ValueError,"rank 2 or 4 is required."
3387          return (arg+transpose(arg))/2
3388        elif isinstance(arg,escript.Data):
3389          if arg.getRank()==2:
3390            if not (arg.getShape()[0]==arg.getShape()[1]):
3391               raise ValueError,"argument must be square."
3392            return arg._symmetric()
3393          elif arg.getRank()==4:
3394            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3395               raise ValueError,"argument must be square."
3396            return arg._symmetric()
3397          else:
3398            raise ValueError,"rank 2 or 4 is required."
3399        elif isinstance(arg,float):
3400          return arg
3401        elif isinstance(arg,int):
3402          return float(arg)
3403        elif isinstance(arg,Symbol):
3404          if arg.getRank()==2:
3405            if not (arg.getShape()[0]==arg.getShape()[1]):
3406               raise ValueError,"argument must be square."
3407          elif arg.getRank()==4:
3408            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3409               raise ValueError,"argument must be square."
3410          else:
3411            raise ValueError,"rank 2 or 4 is required."
3412          return (arg+transpose(arg))/2
3413        else:
3414          raise TypeError,"symmetric: Unknown argument type."
3415    
3416    def nonsymmetric(arg):
3417        """
3418        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3419    
3420        @param arg: square matrix. Must have rank 2 or 4 and be square.
3421        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3422        @return: nonsymmetric part of arg
3423        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3424        """
3425        if isinstance(arg,numarray.NumArray):
3426          if arg.rank==2:
3427            if not (arg.shape[0]==arg.shape[1]):
3428               raise ValueError,"nonsymmetric: argument must be square."
3429          elif arg.rank==4:
3430            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3431               raise ValueError,"nonsymmetric: argument must be square."
3432          else:
3433            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3434          return (arg-transpose(arg))/2
3435        elif isinstance(arg,escript.Data):
3436          if arg.getRank()==2:
3437            if not (arg.getShape()[0]==arg.getShape()[1]):
3438               raise ValueError,"argument must be square."
3439            return arg._nonsymmetric()
3440          elif arg.getRank()==4:
3441            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3442               raise ValueError,"argument must be square."
3443            return arg._nonsymmetric()
3444          else:
3445            raise ValueError,"rank 2 or 4 is required."
3446        elif isinstance(arg,float):
3447          return arg
3448        elif isinstance(arg,int):
3449          return float(arg)
3450        elif isinstance(arg,Symbol):
3451          if arg.getRank()==2:
3452            if not (arg.getShape()[0]==arg.getShape()[1]):
3453               raise ValueError,"nonsymmetric: argument must be square."
3454          elif arg.getRank()==4:
3455            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3456               raise ValueError,"nonsymmetric: argument must be square."
3457          else:
3458            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3459          return (arg-transpose(arg))/2
3460        else:
3461          raise TypeError,"nonsymmetric: Unknown argument type."
3462    
3463    def inverse(arg):
3464        """
3465        returns the inverse of the square matrix arg.
3466    
3467        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3468        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3469        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3470        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3471        @note: for L{escript.Data} objects the dimension is restricted to 3.
3472        """
3473        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3474        if isinstance(arg,numarray.NumArray):
3475          return numarray.linear_algebra.inverse(arg)
3476        elif isinstance(arg,escript.Data):
3477          return escript_inverse(arg)
3478        elif isinstance(arg,float):
3479          return 1./arg
3480        elif isinstance(arg,int):
3481          return 1./float(arg)
3482        elif isinstance(arg,Symbol):
3483          return Inverse_Symbol(arg)
3484        else:
3485          raise TypeError,"inverse: Unknown argument type."
3486    
3487    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3488          "arg is a Data objects!!!"
3489          if not arg.getRank()==2:
3490            raise ValueError,"escript_inverse: argument must have rank 2"
3491          s=arg.getShape()      
3492          if not s[0] == s[1]:
3493            raise ValueError,"escript_inverse: argument must be a square matrix."
3494          out=escript.Data(0.,s,arg.getFunctionSpace())
3495          if s[0]==1:
3496              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3497                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3498              out[0,0]=1./arg[0,0]
3499          elif s[0]==2:
3500              A11=arg[0,0]
3501              A12=arg[0,1]
3502              A21=arg[1,0]
3503              A22=arg[1,1]
3504              D = A11*A22-A12*A21
3505              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3506                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3507              D=1./D
3508              out[0,0]= A22*D
3509              out[1,0]=-A21*D
3510              out[0,1]=-A12*D
3511              out[1,1]= A11*D
3512          elif s[0]==3:
3513              A11=arg[0,0]
3514              A21=arg[1,0]
3515              A31=arg[2,0]
3516              A12=arg[0,1]
3517              A22=arg[1,1]
3518              A32=arg[2,1]
3519              A13=arg[0,2]
3520              A23=arg[1,2]
3521              A33=arg[2,2]
3522              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3523              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3524                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3525              D=1./D
3526              out[0,0]=(A22*A33-A23*A32)*D
3527              out[1,0]=(A31*A23-A21*A33)*D
3528              out[2,0]=(A21*A32-A31*A22)*D
3529              out[0,1]=(A13*A32-A12*A33)*D
3530              out[1,1]=(A11*A33-A31*A13)*D
3531              out[2,1]=(A12*A31-A11*A32)*D
3532              out[0,2]=(A12*A23-A13*A22)*D
3533              out[1,2]=(A13*A21-A11*A23)*D
3534              out[2,2]=(A11*A22-A12*A21)*D
3535          else:
3536             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3537          return out
3538    
3539    class Inverse_Symbol(DependendSymbol):
3540       """
3541       L{Symbol} representing the result of the inverse function
3542       """
3543       def __init__(self,arg):
3544          """
3545          initialization of inverse L{Symbol} with argument arg
3546          @param arg: argument of function
3547          @type arg: L{Symbol}.
3548          """
3549          if not arg.getRank()==2:
3550            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3551          s=arg.getShape()
3552          if not s[0] == s[1]:
3553            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3554          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3555    
3556       def getMyCode(self,argstrs,format="escript"):
3557          """
3558          returns a program code that can be used to evaluate the symbol.
3559    
3560          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3561          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3562          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3563          @type format: C{str}
3564          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3565          @rtype: C{str}
3566          @raise NotImplementedError: if the requested format is not available
3567          """
3568          if format=="escript" or format=="str"  or format=="text":
3569             return "inverse(%s)"%argstrs[0]
3570          else:
3571             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3572    
3573       def substitute(self,argvals):
3574          """
3575          assigns new values to symbols in the definition of the symbol.
3576          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3577    
3578          @param argvals: new values assigned to symbols
3579          @type argvals: C{dict} with keywords of type L{Symbol}.
3580          @return: result of the substitution process. Operations are executed as much as possible.
3581          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3582          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3583          """
3584          if argvals.has_key(self):
3585             arg=argvals[self]
3586             if self.isAppropriateValue(arg):
3587                return arg
3588             else:
3589                raise TypeError,"%s: new value is not appropriate."%str(self)
3590          else:
3591             arg=self.getSubstitutedArguments(argvals)
3592             return inverse(arg[0])
3593    
3594       def diff(self,arg):
3595          """
3596          differential of this object
3597    
3598          @param arg: the derivative is calculated with respect to arg
3599          @type arg: L{escript.Symbol}
3600          @return: derivative with respect to C{arg}
3601          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3602          """
3603          if arg==self:
3604             return identity(self.getShape())
3605          else:
3606             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3607    
3608    def eigenvalues(arg):
3609        """
3610        returns the eigenvalues of the square matrix arg.
3611    
3612        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3613                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3614        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3615        @return: the eigenvalues in increasing order.
3616        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3617        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3618        """
3619        if isinstance(arg,numarray.NumArray):
3620          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3621          out.sort()
3622          return out
3623        elif isinstance(arg,escript.Data):
3624          return arg._eigenvalues()
3625        elif isinstance(arg,Symbol):
3626          if not arg.getRank()==2:
3627            raise ValueError,"eigenvalues: argument must have rank 2"
3628          s=arg.getShape()      
3629          if not s[0] == s[1]:
3630            raise ValueError,"eigenvalues: argument must be a square matrix."
3631          if s[0]==1:
3632              return arg[0]
3633          elif s[0]==2:
3634              arg1=symmetric(arg)
3635              A11=arg1[0,0]
3636              A12=arg1[0,1]
3637              A22=arg1[1,1]
3638              trA=(A11+A22)/2.
3639              A11-=trA
3640              A22-=trA
3641              s=sqrt(A12**2-A11*A22)
3642              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3643          elif s[0]==3:
3644              arg1=symmetric(arg)
3645              A11=arg1[0,0]
3646              A12=arg1[0,1]
3647              A22=arg1[1,1]
3648              A13=arg1[0,2]
3649              A23=arg1[1,2]
3650              A33=arg1[2,2]
3651              trA=(A11+A22+A33)/3.
3652              A11-=trA
3653              A22-=trA
3654              A33-=trA
3655              A13_2=A13**2
3656              A23_2=A23**2
3657              A12_2=A12**2
3658              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3659              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3660              sq_p=sqrt(p/3.)
3661              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
3662              sq_p*=2.
3663              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3664               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3665               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3666              return trA+sq_p*f
3667          else:
3668             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3669        elif isinstance(arg,float):
3670          return arg
3671        elif isinstance(arg,int):
3672          return float(arg)
3673        else:
3674          raise TypeError,"eigenvalues: Unknown argument type."
3675    
3676    def eigenvalues_and_eigenvectors(arg):
3677        """
3678        returns the eigenvalues and eigenvectors of the square matrix arg.
3679    
3680        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3681                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3682        @type arg: L{escript.Data}
3683        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3684                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3685                 the eigenvector coresponding to the i-th eigenvalue.
3686        @rtype: L{tuple} of L{escript.Data}.
3687        @note: The dimension is restricted to 3.
3688        """
3689        if isinstance(arg,numarray.NumArray):
3690          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3691        elif isinstance(arg,escript.Data):
3692          return arg._eigenvalues_and_eigenvectors()
3693        elif isinstance(arg,Symbol):
3694          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3695        elif isinstance(arg,float):
3696          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3697        elif isinstance(arg,int):
3698          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3699        else:
3700          raise TypeError,"eigenvalues: Unknown argument type."
3701  #=======================================================  #=======================================================
3702  #  Binary operations:  #  Binary operations:
3703  #=======================================================  #=======================================================
# Line 2920  class Add_Symbol(DependendSymbol): Line 3741  class Add_Symbol(DependendSymbol):
3741         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3742         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3743         """         """
3744         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3745         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3746         if not sh0==sh1:         if not sh0==sh1:
3747            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3748         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 2936  class Add_Symbol(DependendSymbol): Line 3757  class Add_Symbol(DependendSymbol):
3757        @type format: C{str}        @type format: C{str}
3758        @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.
3759        @rtype: C{str}        @rtype: C{str}
3760        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3761        """        """
3762        if format=="str" or format=="text":        if format=="str" or format=="text":
3763           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 2995  def mult(arg0,arg1): Line 3816  def mult(arg0,arg1):
3816         """         """
3817         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3818         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3819            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3820         else:         else:
3821            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3822                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3019  class Mult_Symbol(DependendSymbol): Line 3840  class Mult_Symbol(DependendSymbol):
3840         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3841         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3842         """         """
3843         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3844         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3845         if not sh0==sh1:         if not sh0==sh1:
3846            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3847         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3035  class Mult_Symbol(DependendSymbol): Line 3856  class Mult_Symbol(DependendSymbol):
3856        @type format: C{str}        @type format: C{str}
3857        @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.
3858        @rtype: C{str}        @rtype: C{str}
3859        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3860        """        """
3861        if format=="str" or format=="text":        if format=="str" or format=="text":
3862           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3095  def quotient(arg0,arg1): Line 3916  def quotient(arg0,arg1):
3916         """         """
3917         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3918         if testForZero(args[0]):         if testForZero(args[0]):
3919            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3920         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3921            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3922               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3124  class Quotient_Symbol(DependendSymbol): Line 3945  class Quotient_Symbol(DependendSymbol):
3945         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3946         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3947         """         """
3948         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3949         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3950         if not sh0==sh1:         if not sh0==sh1:
3951            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3952         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3140  class Quotient_Symbol(DependendSymbol): Line 3961  class Quotient_Symbol(DependendSymbol):
3961        @type format: C{str}        @type format: C{str}
3962        @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.
3963        @rtype: C{str}        @rtype: C{str}
3964        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3965        """        """
3966        if format=="str" or format=="text":        if format=="str" or format=="text":
3967           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3201  def power(arg0,arg1): Line 4022  def power(arg0,arg1):
4022         """         """
4023         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4024         if testForZero(args[0]):         if testForZero(args[0]):
4025            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4026         elif testForZero(args[1]):         elif testForZero(args[1]):
4027            return numarray.ones(args[0],numarray.Float)            return numarray.ones(getShape(args[1]),numarray.Float64)
4028         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4029            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4030         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):
# Line 3226  class Power_Symbol(DependendSymbol): Line 4047  class Power_Symbol(DependendSymbol):
4047         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4048         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4049         """         """
4050         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4051         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4052         if not sh0==sh1:         if not sh0==sh1:
4053            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4054         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3244  class Power_Symbol(DependendSymbol): Line 4065  class Power_Symbol(DependendSymbol):
4065        @type format: C{str}        @type format: C{str}
4066        @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.
4067        @rtype: C{str}        @rtype: C{str}
4068        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4069        """        """
4070        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4071           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3304  def maximum(*args): Line 4125  def maximum(*args):
4125         if out==None:         if out==None:
4126            out=a            out=a
4127         else:         else:
4128            m=whereNegative(out-a)            diff=add(a,-out)
4129            out=m*a+(1.-m)*out            out=add(out,mult(wherePositive(diff),diff))
4130      return out      return out
4131        
4132  def minimum(*arg):  def minimum(*args):
4133      """      """
4134      the minimum over arguments args      the minimum over arguments args
4135    
# Line 3322  def minimum(*arg): Line 4143  def minimum(*arg):
4143         if out==None:         if out==None:
4144            out=a            out=a
4145         else:         else:
4146            m=whereNegative(out-a)            diff=add(a,-out)
4147            out=m*out+(1.-m)*a            out=add(out,mult(whereNegative(diff),diff))
4148      return out      return out
4149    
4150    def clip(arg,minval=None,maxval=None):
4151        """
4152        cuts the values of arg between minval and maxval
4153    
4154        @param arg: argument
4155        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4156        @param minval: lower range. If None no lower range is applied
4157        @type minval: C{float} or C{None}
4158        @param maxval: upper range. If None no upper range is applied
4159        @type maxval: C{float} or C{None}
4160        @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4161                 less then maxval are unchanged.
4162        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4163        @raise ValueError: if minval>maxval
4164        """
4165        if not minval==None and not maxval==None:
4166           if minval>maxval:
4167              raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4168        if minval == None:
4169            tmp=arg
4170        else:
4171            tmp=maximum(minval,arg)
4172        if maxval == None:
4173            return tmp
4174        else:
4175            return minimum(tmp,maxval)
4176    
4177        
4178  def inner(arg0,arg1):  def inner(arg0,arg1):
4179      """      """
# Line 3340  def inner(arg0,arg1): Line 4189  def inner(arg0,arg1):
4189      @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}
4190      @param arg1: second argument      @param arg1: second argument
4191      @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}
4192      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4193      @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
4194      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4195      """      """
4196      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4197      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4198      if not sh0==sh1:      if not sh0==sh1:
4199          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4200      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4201    
4202    def outer(arg0,arg1):
4203        """
4204        the outer product of the two argument:
4205    
4206        out[t,s]=arg0[t]*arg1[s]
4207    
4208        where
4209    
4210            - s runs through arg0.Shape
4211            - t runs through arg1.Shape
4212    
4213        @param arg0: first argument
4214        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4215        @param arg1: second argument
4216        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4217        @return: the outer product of arg0 and arg1 at each data point
4218        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4219        """
4220        return generalTensorProduct(arg0,arg1,axis_offset=0)
4221    
4222  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4223      """      """
4224        see L{matrix_mult}
4225        """
4226        return matrix_mult(arg0,arg1)
4227    
4228    def matrix_mult(arg0,arg1):
4229        """
4230      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4231    
4232      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4233    
4234            or      or
4235    
4236      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4237    
4238      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.
4239    
4240      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4241      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3370  def matrixmult(arg0,arg1): Line 4245  def matrixmult(arg0,arg1):
4245      @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
4246      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4247      """      """
4248      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4249      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4250      if not len(sh0)==2 :      if not len(sh0)==2 :
4251          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4252      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4253          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4254      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4255    
4256  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4257      """      """
4258      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  
4259      """      """
4260      return generalTensorProduct(arg0,arg1,offset=0)      return tensor_mult(arg0,arg1)
4261    
4262    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4263      """      """
4264      the tensor product of the two argument:      the tensor product of the two argument:
   
4265            
4266      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4267    
4268      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4269    
4270                   or      or
4271    
4272      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4273    
# Line 3415  def tensormult(arg0,arg1): Line 4276  def tensormult(arg0,arg1):
4276    
4277      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]
4278                                
4279                   or      or
4280    
4281      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]
4282    
4283                   or      or
4284    
4285      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]
4286    
4287      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  
4288      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.
4289    
4290      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4291      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3433  def tensormult(arg0,arg1): Line 4294  def tensormult(arg0,arg1):
4294      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4295      @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
4296      """      """
4297      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4298      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4299      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4300         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4301      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):
4302         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4303      else:      else:
4304          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4305    
4306  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4307      """      """
4308      generalized tensor product      generalized tensor product
4309    
4310      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4311    
4312      where s runs through arg0.Shape[:arg0.Rank-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[offset:]  
4313    
4314      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]
4315      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4316            - t runs through arg1.Shape[axis_offset:]
4317    
4318      @param arg0: first argument      @param arg0: first argument
4319      @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}
4320      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4321      @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}
4322      @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.
4323      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
# Line 3467  def generalTensorProduct(arg0,arg1,offse Line 4327  def generalTensorProduct(arg0,arg1,offse
4327      # 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
4328      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4329         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4330             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4331         else:         else:
4332             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4333                 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)
4334             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4335             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4336             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4337             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4338             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4339             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4340             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4341             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4342             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4343             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4344             for i0 in range(d0):             for i0 in range(d0):
4345                      for i1 in range(d1):                      for i1 in range(d1):
4346                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4347             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4348             return out             return out
4349      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4350         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4351             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4352         else:         else:
4353             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)
4354      else:            else:      
4355         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4356                                    
4357  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4358     """     """
4359     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4360     """     """
4361     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4362         """         """
4363         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4364    
4365         @param arg0: numerator         @param arg0: first argument
4366         @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}.
4367         @param arg1: denominator         @param arg1: second argument
4368         @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}.
4369         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4370         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4371         """         """
4372         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4373         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4374         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4375         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4376         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4377         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4378         if not sh01==sh10:         if not sh01==sh10:
4379             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)
4380         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])
4381    
4382     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4383        """        """
# Line 3529  class GeneralTensorProduct_Symbol(Depend Line 4389  class GeneralTensorProduct_Symbol(Depend
4389        @type format: C{str}        @type format: C{str}
4390        @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.
4391        @rtype: C{str}        @rtype: C{str}
4392        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4393        """        """
4394        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4395           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])
4396        else:        else:
4397           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)
4398    
# Line 3557  class GeneralTensorProduct_Symbol(Depend Line 4417  class GeneralTensorProduct_Symbol(Depend
4417           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4418           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4419    
4420  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4421      "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!!!"
4422      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4423      shape0=arg0.getShape()[:arg0.getRank()-offset]  
4424      shape01=arg0.getShape()[arg0.getRank()-offset:]  def transposed_matrix_mult(arg0,arg1):
4425      shape10=arg1.getShape()[:offset]      """
4426      shape1=arg1.getShape()[offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4427      if not shape01==shape10:  
4428          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]
4429    
4430      # create return value:      or
4431      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())  
4432      #      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4433      s0=[[]]  
4434      for k in shape0:      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4435            s=[]  
4436            for j in s0:      The first dimension of arg0 and arg1 must match.
4437                  for i in range(k): s.append(j+[slice(i,i)])  
4438            s0=s      @param arg0: first argument of rank 2
4439      s1=[[]]      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4440      for k in shape1:      @param arg1: second argument of at least rank 1
4441            s=[]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4442            for j in s1:      @return: the product of the transposed of arg0 and arg1 at each data point
4443                  for i in range(k): s.append(j+[slice(i,i)])      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4444            s1=s      @raise ValueError: if the shapes of the arguments are not appropriate
4445      s01=[[]]      """
4446      for k in shape01:      sh0=getShape(arg0)
4447            s=[]      sh1=getShape(arg1)
4448            for j in s01:      if not len(sh0)==2 :
4449                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"first argument must have rank 2"
4450            s01=s      if not len(sh1)==2 and not len(sh1)==1:
4451            raise ValueError,"second argument must have rank 1 or 2"
4452      for i0 in s0:      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4453         for i1 in s1:  
4454           s=escript.Scalar(0.,arg0.getFunctionSpace())  def transposed_tensor_mult(arg0,arg1):
4455           for i01 in s01:      """
4456              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      the tensor product of the transposed of the first and the second argument
4457           out.__setitem__(tuple(i0+i1),s)      
4458      return out      for arg0 of rank 2 this is
4459    
4460        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4461    
4462        or
4463    
4464        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4465    
4466      
4467        and for arg0 of rank 4 this is
4468    
4469        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4470                  
4471        or
4472    
4473        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4474    
4475        or
4476    
4477        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4478    
4479        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4480        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4481    
4482        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4483    
4484        @param arg0: first argument of rank 2 or 4
4485        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4486        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4487        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4488        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4489        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4490        """
4491        sh0=getShape(arg0)
4492        sh1=getShape(arg1)
4493        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4494           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4495        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4496           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4497        else:
4498            raise ValueError,"first argument must have rank 2 or 4"
4499    
4500    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4501        """
4502        generalized tensor product of transposed of arg0 and arg1:
4503    
4504        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4505    
4506        where
4507    
4508            - s runs through arg0.Shape[axis_offset:]
4509            - r runs trough arg0.Shape[:axis_offset]
4510            - t runs through arg1.Shape[axis_offset:]
4511    
4512        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4513        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4514    
4515        @param arg0: first argument
4516        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4517        @param arg1: second argument
4518        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4519        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4520        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4521        """
4522        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4523        arg0,arg1=matchType(arg0,arg1)
4524        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4525        if isinstance(arg0,numarray.NumArray):
4526           if isinstance(arg1,Symbol):
4527               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4528           else:
4529               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4530                   raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4531               arg0_c=arg0.copy()
4532               arg1_c=arg1.copy()
4533               sh0,sh1=arg0.shape,arg1.shape
4534               d0,d1,d01=1,1,1
4535               for i in sh0[axis_offset:]: d0*=i
4536               for i in sh1[axis_offset:]: d1*=i
4537               for i in sh0[:axis_offset]: d01*=i
4538               arg0_c.resize((d01,d0))
4539               arg1_c.resize((d01,d1))
4540               out=numarray.zeros((d0,d1),numarray.Float64)
4541               for i0 in range(d0):
4542                        for i1 in range(d1):
4543                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4544               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4545               return out
4546        elif isinstance(arg0,escript.Data):
4547           if isinstance(arg1,Symbol):
4548               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4549           else:
4550               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4551        else:      
4552           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4553                    
4554    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4555       """
4556       Symbol representing the general tensor product of the transposed of arg0 and arg1
4557       """
4558       def __init__(self,arg0,arg1,axis_offset=0):
4559           """
4560           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4561    
4562           @param arg0: first argument
4563           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4564           @param arg1: second argument
4565           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4566           @raise ValueError: inconsistent dimensions of arguments.
4567           @note: if both arguments have a spatial dimension, they must equal.
4568           """
4569           sh_arg0=getShape(arg0)
4570           sh_arg1=getShape(arg1)
4571           sh01=sh_arg0[:axis_offset]
4572           sh10=sh_arg1[:axis_offset]
4573           sh0=sh_arg0[axis_offset:]
4574           sh1=sh_arg1[axis_offset:]
4575           if not sh01==sh10:
4576               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)
4577           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4578    
4579       def getMyCode(self,argstrs,format="escript"):
4580          """
4581          returns a program code that can be used to evaluate the symbol.
4582    
4583          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4584          @type argstrs: C{list} of length 2 of C{str}.
4585          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4586          @type format: C{str}
4587          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4588          @rtype: C{str}
4589          @raise NotImplementedError: if the requested format is not available
4590          """
4591          if format=="escript" or format=="str" or format=="text":
4592             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4593          else:
4594             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4595    
4596       def substitute(self,argvals):
4597          """
4598          assigns new values to symbols in the definition of the symbol.
4599          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4600    
4601          @param argvals: new values assigned to symbols
4602          @type argvals: C{dict} with keywords of type L{Symbol}.
4603          @return: result of the substitution process. Operations are executed as much as possible.
4604          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4605          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4606          """
4607          if argvals.has_key(self):
4608             arg=argvals[self]
4609             if self.isAppropriateValue(arg):
4610                return arg
4611             else:
4612                raise TypeError,"%s: new value is not appropriate."%str(self)
4613          else:
4614             args=self.getSubstitutedArguments(argvals)
4615             return generalTransposedTensorProduct(args[0],args[1],args[2])
4616    
4617    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4618        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4619        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4620    
4621    def matrix_transposed_mult(arg0,arg1):
4622        """
4623        matrix-transposed(matrix) product of the two argument:
4624    
4625        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4626    
4627        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4628    
4629        The last dimensions of arg0 and arg1 must match.
4630    
4631        @param arg0: first argument of rank 2
4632        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4633        @param arg1: second argument of rank 2
4634        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4635        @return: the product of arg0 and the transposed of arg1 at each data point
4636        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4637        @raise ValueError: if the shapes of the arguments are not appropriate
4638        """
4639        sh0=getShape(arg0)
4640        sh1=getShape(arg1)
4641        if not len(sh0)==2 :
4642            raise ValueError,"first argument must have rank 2"
4643        if not len(sh1)==2 and not len(sh1)==1:
4644            raise ValueError,"second argument must have rank 1 or 2"
4645        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4646    
4647    def tensor_transposed_mult(arg0,arg1):
4648        """
4649        the tensor product of the first and the transpose of the second argument
4650        
4651        for arg0 of rank 2 this is
4652    
4653        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4654    
4655        and for arg0 of rank 4 this is
4656    
4657        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4658                  
4659        or
4660    
4661        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4662    
4663        In the first case the the second dimension of arg0 and arg1 must match and  
4664        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4665    
4666        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4667    
4668        @param arg0: first argument of rank 2 or 4
4669        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4670        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4671        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4672        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4673        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4674        """
4675        sh0=getShape(arg0)
4676        sh1=getShape(arg1)
4677        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4678           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4679        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4680           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4681        else:
4682            raise ValueError,"first argument must have rank 2 or 4"
4683    
4684    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4685        """
4686        generalized tensor product of transposed of arg0 and arg1:
4687    
4688        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4689    
4690        where
4691    
4692            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4693            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4694            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4695    
4696        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4697        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4698    
4699        @param arg0: first argument
4700        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4701        @param arg1: second argument
4702        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4703        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4704        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4705        """
4706        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4707        arg0,arg1=matchType(arg0,arg1)
4708        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4709        if isinstance(arg0,numarray.NumArray):
4710           if isinstance(arg1,Symbol):
4711               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4712           else:
4713               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4714                   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)
4715               arg0_c=arg0.copy()
4716               arg1_c=arg1.copy()
4717               sh0,sh1=arg0.shape,arg1.shape
4718               d0,d1,d01=1,1,1
4719               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4720               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4721               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4722               arg0_c.resize((d0,d01))
4723               arg1_c.resize((d1,d01))
4724               out=numarray.zeros((d0,d1),numarray.Float64)
4725               for i0 in range(d0):
4726                        for i1 in range(d1):
4727                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4728               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4729               return out
4730        elif isinstance(arg0,escript.Data):
4731           if isinstance(arg1,Symbol):
4732               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4733           else:
4734               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4735        else:      
4736           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4737                    
4738    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4739       """
4740       Symbol representing the general tensor product of arg0 and the transpose of arg1
4741       """
4742       def __init__(self,arg0,arg1,axis_offset=0):
4743           """
4744           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4745    
4746           @param arg0: first argument
4747           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4748           @param arg1: second argument
4749           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4750           @raise ValueError: inconsistent dimensions of arguments.
4751           @note: if both arguments have a spatial dimension, they must equal.
4752           """
4753           sh_arg0=getShape(arg0)
4754           sh_arg1=getShape(arg1)
4755           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4756           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4757           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4758           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4759           if not sh01==sh10:
4760               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)
4761           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4762    
4763       def getMyCode(self,argstrs,format="escript"):
4764          """
4765          returns a program code that can be used to evaluate the symbol.
4766    
4767          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4768          @type argstrs: C{list} of length 2 of C{str}.
4769          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4770          @type format: C{str}
4771          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4772          @rtype: C{str}
4773          @raise NotImplementedError: if the requested format is not available
4774          """
4775          if format=="escript" or format=="str" or format=="text":
4776             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4777          else:
4778             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4779    
4780       def substitute(self,argvals):
4781          """
4782          assigns new values to symbols in the definition of the symbol.
4783          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4784    
4785          @param argvals: new values assigned to symbols
4786          @type argvals: C{dict} with keywords of type L{Symbol}.
4787          @return: result of the substitution process. Operations are executed as much as possible.
4788          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4789          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4790          """
4791          if argvals.has_key(self):
4792             arg=argvals[self]
4793             if self.isAppropriateValue(arg):
4794                return arg
4795             else:
4796                raise TypeError,"%s: new value is not appropriate."%str(self)
4797          else:
4798             args=self.getSubstitutedArguments(argvals)
4799             return generalTensorTransposedProduct(args[0],args[1],args[2])
4800    
4801    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4802        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4803        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4804    
4805  #=========================================================  #=========================================================
4806  #   some little helpers  #  functions dealing with spatial dependency
4807  #=========================================================  #=========================================================
4808  def grad(arg,where=None):  def grad(arg,where=None):
4809      """      """
4810      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4811    
4812        If C{g} is the returned object, then
4813    
4814      @param arg:   Data object representing the function which gradient        - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4815                    to be calculated.        - if C{arg} is rank 1 C{g[i,s]} is the derivative of C{arg[i]} with respect to the C{s}-th spatial dimension.
4816          - 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.
4817          - 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.
4818    
4819        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4820        @type arg: L{escript.Data} or L{Symbol}
4821      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4822                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4823        @type where: C{None} or L{escript.FunctionSpace}
4824        @return: gradient of arg.
4825        @rtype: L{escript.Data} or L{Symbol}
4826      """      """
4827      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4828         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3617  def grad(arg,where=None): Line 4832  def grad(arg,where=None):
4832         else:         else:
4833            return arg._grad(where)            return arg._grad(where)
4834      else:      else:
4835        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4836    
4837    class Grad_Symbol(DependendSymbol):
4838       """
4839       L{Symbol} representing the result of the gradient operator
4840       """
4841       def __init__(self,arg,where=None):
4842          """
4843          initialization of gradient L{Symbol} with argument arg
4844          @param arg: argument of function
4845          @type arg: L{Symbol}.
4846          @param where: FunctionSpace in which the gradient will be calculated.
4847                      If not present or C{None} an appropriate default is used.
4848          @type where: C{None} or L{escript.FunctionSpace}
4849          """
4850          d=arg.getDim()
4851          if d==None:
4852             raise ValueError,"argument must have a spatial dimension"
4853          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4854    
4855       def getMyCode(self,argstrs,format="escript"):
4856          """
4857          returns a program code that can be used to evaluate the symbol.
4858    
4859          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4860          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4861          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4862          @type format: C{str}
4863          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4864          @rtype: C{str}
4865          @raise NotImplementedError: if the requested format is not available
4866          """
4867          if format=="escript" or format=="str"  or format=="text":
4868             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4869          else:
4870             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4871    
4872       def substitute(self,argvals):
4873          """
4874          assigns new values to symbols in the definition of the symbol.
4875          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4876    
4877          @param argvals: new values assigned to symbols
4878          @type argvals: C{dict} with keywords of type L{Symbol}.
4879          @return: result of the substitution process. Operations are executed as much as possible.
4880          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4881          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4882          """
4883          if argvals.has_key(self):
4884             arg=argvals[self]
4885             if self.isAppropriateValue(arg):
4886                return arg
4887             else:
4888                raise TypeError,"%s: new value is not appropriate."%str(self)
4889          else:
4890             arg=self.getSubstitutedArguments(argvals)
4891             return grad(arg[0],where=arg[1])
4892    
4893       def diff(self,arg):
4894          """
4895          differential of this object
4896    
4897          @param arg: the derivative is calculated with respect to arg
4898          @type arg: L{escript.Symbol}
4899          @return: derivative with respect to C{arg}
4900          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4901          """
4902          if arg==self:
4903             return identity(self.getShape())
4904          else:
4905             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4906    
4907  def integrate(arg,where=None):  def integrate(arg,where=None):
4908      """      """
4909      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}
4910      its domain.      before integration.
4911    
4912      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4913        @type arg: L{escript.Data} or L{Symbol}
4914      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4915                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4916        @type where: C{None} or L{escript.FunctionSpace}
4917        @return: integral of arg.
4918        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4919      """      """
4920      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):  
4921         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4922      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4923         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
# Line 3652  def integrate(arg,where=None): Line 4926  def integrate(arg,where=None):
4926         else:         else:
4927            return arg._integrate()            return arg._integrate()
4928      else:      else:
4929        raise TypeError,"integrate: Unknown argument type."         arg2=escript.Data(arg,where)
4930           if arg2.getRank()==0:
4931              return arg2._integrate()[0]
4932           else:
4933              return arg2._integrate()
4934    
4935    class Integrate_Symbol(DependendSymbol):
4936       """
4937       L{Symbol} representing the result of the spatial integration operator
4938       """
4939       def __init__(self,arg,where=None):
4940          """
4941          initialization of integration L{Symbol} with argument arg
4942          @param arg: argument of the integration
4943          @type arg: L{Symbol}.
4944          @param where: FunctionSpace in which the integration will be calculated.
4945                      If not present or C{None} an appropriate default is used.
4946          @type where: C{None} or L{escript.FunctionSpace}
4947          """
4948          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4949    
4950       def getMyCode(self,argstrs,format="escript"):
4951          """
4952          returns a program code that can be used to evaluate the symbol.
4953    
4954          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4955          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4956          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4957          @type format: C{str}
4958          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4959          @rtype: C{str}
4960          @raise NotImplementedError: if the requested format is not available
4961          """
4962          if format=="escript" or format=="str"  or format=="text":
4963             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4964          else:
4965             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4966    
4967       def substitute(self,argvals):
4968          """
4969          assigns new values to symbols in the definition of the symbol.
4970          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4971    
4972          @param argvals: new values assigned to symbols
4973          @type argvals: C{dict} with keywords of type L{Symbol}.
4974          @return: result of the substitution process. Operations are executed as much as possible.
4975          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4976          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4977          """
4978          if argvals.has_key(self):
4979             arg=argvals[self]
4980             if self.isAppropriateValue(arg):
4981                return arg
4982             else:
4983                raise TypeError,"%s: new value is not appropriate."%str(self)
4984          else:
4985             arg=self.getSubstitutedArguments(argvals)
4986             return integrate(arg[0],where=arg[1])
4987    
4988       def diff(self,arg):
4989          """
4990          differential of this object
4991    
4992          @param arg: the derivative is calculated with respect to arg
4993          @type arg: L{escript.Symbol}
4994          @return: derivative with respect to C{arg}
4995          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4996          """
4997          if arg==self:
4998             return identity(self.getShape())
4999          else:
5000             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5001    
5002    
5003  def interpolate(arg,where):  def interpolate(arg,where):
5004      """      """
5005      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
5006    
5007      @param arg:    interpolant      @param arg: interpolant
5008      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
5009        @param where: FunctionSpace to be interpolated to
5010        @type where: L{escript.FunctionSpace}
5011        @return: interpolated argument
5012        @rtype: C{escript.Data} or L{Symbol}
5013      """      """
5014      if testForZero(arg):      if isinstance(arg,Symbol):
5015        return 0         return Interpolate_Symbol(arg,where)
     elif isinstance(arg,Symbol):  
        return Interpolated_Symbol(arg,where)  
5016      else:      else:
5017         return escript.Data(arg,where)         return escript.Data(arg,where)
5018    
5019  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
5020      """     """
5021      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
5022       """
5023       def __init__(self,arg,where):
5024          """
5025          initialization of interpolation L{Symbol} with argument arg
5026          @param arg: argument of the interpolation
5027          @type arg: L{Symbol}.
5028          @param where: FunctionSpace into which the argument is interpolated.
5029          @type where: L{escript.FunctionSpace}
5030          """
5031          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
5032    
5033      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
5034                    be calculated.        """
5035      @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)  
5036    
5037  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
5038      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
5039      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.
5040          @type format: C{str}
5041          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
5042          @rtype: C{str}
5043          @raise NotImplementedError: if the requested format is not available
5044          """
5045          if format=="escript" or format=="str"  or format=="text":
5046             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
5047          else:
5048             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
5049    
5050      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
5051                    to be calculated.        """
5052      """        assigns new values to symbols in the definition of the symbol.
5053      d=arg.getDomain()        The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
     return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())  
5054    
5055  #=============================        @param argvals: new values assigned to symbols
5056  #        @type argvals: C{dict} with keywords of type L{Symbol}.
5057  # 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.
5058  # 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
5059  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
5060          """
5061          if argvals.has_key(self):
5062             arg=argvals[self]
5063             if self.isAppropriateValue(arg):
5064                return arg
5065             else:
5066                raise TypeError,"%s: new value is not appropriate."%str(self)
5067          else:
5068             arg=self.getSubstitutedArguments(argvals)
5069             return interpolate(arg[0],where=arg[1])
5070    
5071       def diff(self,arg):
5072          """
5073          differential of this object
5074    
5075          @param arg: the derivative is calculated with respect to arg
5076          @type arg: L{escript.Symbol}
5077          @return: derivative with respect to C{arg}
5078          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
5079          """
5080          if arg==self:
5081             return identity(self.getShape())
5082          else:
5083             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5084    
 # functions involving the underlying Domain:  
5085    
5086  def transpose(arg,axis=None):  def div(arg,where=None):
5087      """      """
5088      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
5089    
5090      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
5091        @type arg: L{escript.Data} or L{Symbol}
5092        @param where: FunctionSpace in which the divergence will be calculated.
5093                      If not present or C{None} an appropriate default is used.
5094        @type where: C{None} or L{escript.FunctionSpace}
5095        @return: divergence of arg.
5096        @rtype: L{escript.Data} or L{Symbol}
5097      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
5098      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5099         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
5100      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
5101         # 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)  
5102      else:      else:
5103         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
5104        if not arg.getShape()==(dim,):
5105          raise ValueError,"div: expected shape is (%s,)"%dim
5106        return trace(grad(arg,where))
5107    
5108  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
5109      """      """
5110      Return      returns the jump of arg across the continuity of the domain
5111    
5112      @param arg:      @param arg: argument
5113        @type arg: L{escript.Data} or L{Symbol}
5114        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5115                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5116        @type domain: C{None} or L{escript.Domain}
5117        @return: jump of arg
5118        @rtype: L{escript.Data} or L{Symbol}
5119      """      """
5120      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
5121         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)  
5122    
5123    def L2(arg):
5124        """
5125        returns the L2 norm of arg at where
5126        
5127        @param arg: function which L2 to be calculated.
5128        @type arg: L{escript.Data} or L{Symbol}
5129        @return: L2 norm of arg.
5130        @rtype: L{float} or L{Symbol}
5131        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5132        """
5133        return sqrt(integrate(inner(arg,arg)))
5134    #=============================
5135    #
5136    
5137  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5138      """      """
5139      resorts the component of arg according to index      resorts the component of arg according to index
5140    
5141      """      """
5142      pass      raise NotImplementedError
5143  #  #
5144  # $Log: util.py,v $  # $Log: util.py,v $
5145  # Revision 1.14.2.16  2005/10/19 06:09:57  gross  # Revision 1.14.2.16  2005/10/19 06:09:57  gross

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

  ViewVC Help
Powered by ViewVC 1.1.26