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

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

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

trunk/escript/py_src/util.py revision 341 by gross, Mon Dec 12 05:26:10 2005 UTC temp/escript/py_src/util.py revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC
# Line 1  Line 1 
1    #
2  # $Id$  # $Id$
3  #  #
4  #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  #######################################################
5    #
6    #           Copyright 2003-2007 by ACceSS MNRF
7    #       Copyright 2007 by University of Queensland
8  #  #
9  #   This software is the property of ACcESS.  No part of this code  #                http://esscc.uq.edu.au
10  #   may be copied in any form or by any means without the expressed written  #        Primary Business: Queensland, Australia
11  #   consent of ACcESS.  Copying, use or modification of this software  #  Licensed under the Open Software License version 3.0
12  #   by any unauthorised person is illegal unless that  #     http://www.opensource.org/licenses/osl-3.0.php
13  #   person has a software license agreement with ACcESS.  #
14    #######################################################
15  #  #
16    
17  """  """
18  Utility functions for escript  Utility functions for escript
19    
 @remark:  This module is under construction and is still tested!!!  
   
20  @var __author__: name of author  @var __author__: name of author
21  @var __licence__: licence agreement  @var __copyright__: copyrights
22    @var __license__: licence agreement
23  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
24  @var __version__: version  @var __version__: version
25  @var __date__: date of the version  @var __date__: date of the version
26    @var EPSILON: smallest positive value with 1.<1.+EPSILON
27  """  """
28                                                                                                                                                                                                                                                                                                                                                                                                            
29  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
30  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
31                        http://www.access.edu.au
32                    Primary Business: Queensland, Australia"""
33    __license__="""Licensed under the Open Software License version 3.0
34                 http://www.opensource.org/licenses/osl-3.0.php"""
35  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
36  __version__="$Revision: 329 $"  __version__="$Revision$"
37  __date__="$Date$"  __date__="$Date$"
38    
39    
# Line 32  import math Line 41  import math
41  import numarray  import numarray
42  import escript  import escript
43  import os  import os
44    from esys.escript import C_GeneralTensorProduct
45  # 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  
46    
47  #=========================================================  #=========================================================
48  #   some helpers:  #   some helpers:
49  #=========================================================  #=========================================================
50    def getEpsilon():
51         #     ------------------------------------------------------------------
52         #     Compute EPSILON, the machine precision.  The call to daxpp is
53         #     inTENded to fool compilers that use extra-length registers.
54         #     31 Map 1999: Hardwire EPSILON so the debugger can step thru easily.
55         #     ------------------------------------------------------------------
56         eps    = 2.**(-12)
57         p=2.
58         while p>1.:
59                eps/=2.
60                p=1.+eps
61         return eps*2.
62    
63    EPSILON=getEpsilon()
64    
65    def getTagNames(domain):
66        """
67        returns a list of the tag names used by the domain
68    
69        
70        @param domain: a domain object
71        @type domain: L{escript.Domain}
72        @return: a list of the tag name used by the domain.
73        @rtype: C{list} of C{str}
74        """
75        return [n.strip() for n in domain.showTagNames().split(",") ]
76    
77    def insertTagNames(domain,**kwargs):
78        """
79        inserts tag names into the domain
80    
81        @param domain: a domain object
82        @type domain: C{escript.Domain}
83        @keyword <tag name>: tag key assigned to <tag name>
84        @type <tag name>: C{int}
85        """
86        for  k in kwargs:
87             domain.setTagMap(k,kwargs[k])
88    
89    def insertTaggedValues(target,**kwargs):
90        """
91        inserts tagged values into the tagged using tag names
92    
93        @param target: data to be filled by tagged values
94        @type target: L{escript.Data}
95        @keyword <tag name>: value to be used for <tag name>
96        @type <tag name>: C{float} or {numarray.NumArray}
97        @return: C{target}
98        @rtype: L{escript.Data}
99        """
100        for k in kwargs:
101            target.setTaggedValue(k,kwargs[k])
102        return target
103    
104  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
105      """      """
106      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.
107    
108      Example:      Example::
109    
110         tmp=Scalar(..)         tmp=Scalar(..)
111         v=Vector(..)         v=Vector(..)
112         saveVTK("solution.xml",temperature=tmp,velovity=v)         saveVTK("solution.xml",temperature=tmp,velocity=v)
113    
114      tmp and v are written into "solution.xml" where tmp is named "temperature" and v is named "velovity"      tmp and v are written into "solution.xml" where tmp is named "temperature" and v is named "velocity"
115    
116      @param filename: file name of the output file      @param filename: file name of the output file
117      @type filename: C{str}      @type filename: C{str}
# Line 84  def saveVTK(filename,domain=None,**data) Line 121  def saveVTK(filename,domain=None,**data)
121      @type <name>: L{Data} object.      @type <name>: L{Data} object.
122      @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.
123      """      """
124      if domain==None:      new_data={}
125         for i in data.keys():      for n,d in data.items():
126            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
127                fs=d.getFunctionSpace()
128                domain2=fs.getDomain()
129                if fs == escript.Solution(domain2):
130                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
131                elif fs == escript.ReducedSolution(domain2):
132                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
133                else:
134                   new_data[n]=d
135                if domain==None: domain=domain2
136      if domain==None:      if domain==None:
137          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
138      else:      domain.saveVTK(filename,new_data)
         domain.saveVTK(filename,data)  
139    
140  def saveDX(filename,domain=None,**data):  def saveDX(filename,domain=None,**data):
141      """      """
142      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.
143    
144      Example:      Example::
145    
146         tmp=Scalar(..)         tmp=Scalar(..)
147         v=Vector(..)         v=Vector(..)
148         saveDX("solution.dx",temperature=tmp,velovity=v)         saveDX("solution.dx",temperature=tmp,velocity=v)
149    
150      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity".      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velocity".
151    
152      @param filename: file name of the output file      @param filename: file name of the output file
153      @type filename: C{str}      @type filename: C{str}
# Line 112  def saveDX(filename,domain=None,**data): Line 157  def saveDX(filename,domain=None,**data):
157      @type <name>: L{Data} object.      @type <name>: L{Data} object.
158      @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.
159      """      """
160      if domain==None:      new_data={}
161         for i in data.keys():      for n,d in data.items():
162            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
163                fs=d.getFunctionSpace()
164                domain2=fs.getDomain()
165                if fs == escript.Solution(domain2):
166                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
167                elif fs == escript.ReducedSolution(domain2):
168                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
169                elif fs == escript.ContinuousFunction(domain2):
170                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
171                else:
172                   new_data[n]=d
173                if domain==None: domain=domain2
174      if domain==None:      if domain==None:
175          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
176      else:      domain.saveDX(filename,new_data)
         domain.saveDX(filename,data)  
177    
178  def kronecker(d=3):  def kronecker(d=3):
179     """     """
180     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
181    
182     @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
183     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
184     @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
185     @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}  
186     """     """
187     return identityTensor(d)     return identityTensor(d)
188    
# Line 143  def identity(shape=()): Line 197  def identity(shape=()):
197     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
198     """     """
199     if len(shape)>0:     if len(shape)>0:
200        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
201        if len(shape)==1:        if len(shape)==1:
202            for i0 in range(shape[0]):            for i0 in range(shape[0]):
203               out[i0,i0]=1.               out[i0,i0]=1.
   
204        elif len(shape)==2:        elif len(shape)==2:
205            for i0 in range(shape[0]):            for i0 in range(shape[0]):
206               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 163  def identityTensor(d=3): Line 216  def identityTensor(d=3):
216     return the dxd identity matrix     return the dxd identity matrix
217    
218     @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
219     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
220     @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
221     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
222     """     """
223     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
224        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
225     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
226           return identity((d.getDim(),))
227       else:
228           return identity((d,))
229    
230  def identityTensor4(d=3):  def identityTensor4(d=3):
231     """     """
# Line 178  def identityTensor4(d=3): Line 234  def identityTensor4(d=3):
234     @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
235     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
236     @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
237     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
238     """     """
239     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
240        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
241     return identity((d,d))     elif isinstance(d,escript.Domain):
242           return identity((d.getDim(),d.getDim()))
243       else:
244           return identity((d,d))
245    
246  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
247     """     """
# Line 191  def unitVector(i=0,d=3): Line 250  def unitVector(i=0,d=3):
250     @param i: index     @param i: index
251     @type i: C{int}     @type i: C{int}
252     @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
253     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
254     @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
255     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
256     """     """
257     return kronecker(d)[i]     return kronecker(d)[i]
258    
# Line 249  def inf(arg): Line 308  def inf(arg):
308    
309      @param arg: argument      @param arg: argument
310      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
311      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
312      @rtype: C{float}      @rtype: C{float}
313      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
314      """      """
# Line 268  def inf(arg): Line 327  def inf(arg):
327  #=========================================================================  #=========================================================================
328  #   some little helpers  #   some little helpers
329  #=========================================================================  #=========================================================================
330  def pokeShape(arg):  def getRank(arg):
331        """
332        identifies the rank of its argument
333    
334        @param arg: a given object
335        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
336        @return: the rank of the argument
337        @rtype: C{int}
338        @raise TypeError: if type of arg cannot be processed
339        """
340    
341        if isinstance(arg,numarray.NumArray):
342            return arg.rank
343        elif isinstance(arg,escript.Data):
344            return arg.getRank()
345        elif isinstance(arg,float):
346            return 0
347        elif isinstance(arg,int):
348            return 0
349        elif isinstance(arg,Symbol):
350            return arg.getRank()
351        else:
352          raise TypeError,"getShape: cannot identify shape"
353    def getShape(arg):
354      """      """
355      identifies the shape of its argument      identifies the shape of its argument
356    
# Line 290  def pokeShape(arg): Line 372  def pokeShape(arg):
372      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
373          return arg.getShape()          return arg.getShape()
374      else:      else:
375        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
376    
377  def pokeDim(arg):  def pokeDim(arg):
378      """      """
# Line 313  def commonShape(arg0,arg1): Line 395  def commonShape(arg0,arg1):
395      """      """
396      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.
397    
398      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
399      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
400      @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.
401      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
402      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
403      """      """
404      sh0=pokeShape(arg0)      sh0=getShape(arg0)
405      sh1=pokeShape(arg1)      sh1=getShape(arg1)
406      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
407         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
408               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 420  def commonDim(*args):
420      """      """
421      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.
422    
423      @param *args: given objects      @param args: given objects
424      @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
425               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
426      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 360  def testForZero(arg): Line 442  def testForZero(arg):
442    
443      @param arg: a given object      @param arg: a given object
444      @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}
445      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
446      @rtype : C{bool}      @rtype: C{bool}
447      """      """
448      try:      if isinstance(arg,numarray.NumArray):
449           return not Lsup(arg)>0.
450        elif isinstance(arg,escript.Data):
451           return False
452        elif isinstance(arg,float):
453           return not Lsup(arg)>0.
454        elif isinstance(arg,int):
455         return not Lsup(arg)>0.         return not Lsup(arg)>0.
456      except TypeError:      elif isinstance(arg,Symbol):
457           return False
458        else:
459         return False         return False
460    
461  def matchType(arg0=0.,arg1=0.):  def matchType(arg0=0.,arg1=0.):
# Line 386  def matchType(arg0=0.,arg1=0.): Line 476  def matchType(arg0=0.,arg1=0.):
476         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
477            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
478         elif isinstance(arg1,float):         elif isinstance(arg1,float):
479            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
480         elif isinstance(arg1,int):         elif isinstance(arg1,int):
481            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
482         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
483            pass            pass
484         else:         else:
# Line 412  def matchType(arg0=0.,arg1=0.): Line 502  def matchType(arg0=0.,arg1=0.):
502         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
503            pass            pass
504         elif isinstance(arg1,float):         elif isinstance(arg1,float):
505            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
506         elif isinstance(arg1,int):         elif isinstance(arg1,int):
507            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
508         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
509            pass            pass
510         else:         else:
511            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
512      elif isinstance(arg0,float):      elif isinstance(arg0,float):
513         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
514            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
515         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
516            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
517         elif isinstance(arg1,float):         elif isinstance(arg1,float):
518            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
519            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
520         elif isinstance(arg1,int):         elif isinstance(arg1,int):
521            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
522            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
523         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
524            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
525         else:         else:
526            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
527      elif isinstance(arg0,int):      elif isinstance(arg0,int):
528         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
529            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
530         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
531            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
532         elif isinstance(arg1,float):         elif isinstance(arg1,float):
533            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
534            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
535         elif isinstance(arg1,int):         elif isinstance(arg1,int):
536            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
537            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
538         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
539            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
540         else:         else:
541            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
542      else:      else:
# Line 456  def matchType(arg0=0.,arg1=0.): Line 546  def matchType(arg0=0.,arg1=0.):
546    
547  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
548      """      """
549            return representations of arg0 amd arg1 which ahve the same shape
   
     If shape is not given the shape "largest" shape of args is used.  
550    
551      @param args: a given ob      @param arg0: a given object
552      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
553      @return: True if the argument is identical to zero.      @param arg1: a given object
554      @rtype: C{list} of C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
555        @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
556        @rtype: C{tuple}
557      """      """
558      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
559      sh0=pokeShape(arg0)      sh0=getShape(arg0)
560      sh1=pokeShape(arg1)      sh1=getShape(arg1)
561      if len(sh0)<len(sh):      if len(sh0)<len(sh):
562         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
563      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
564         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
565      else:      else:
566         return arg0,arg1         return arg0,arg1
567  #=========================================================  #=========================================================
# Line 491  class Symbol(object): Line 581  class Symbol(object):
581         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
582         symbols or any other object.         symbols or any other object.
583    
584         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
585         @type arg: C{list}         @type args: C{list}
586         @param shape: the shape         @param shape: the shape
587         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
588         @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 625  class Symbol(object):
625         """         """
626         the shape of the symbol.         the shape of the symbol.
627    
628         @return : the shape of the symbol.         @return: the shape of the symbol.
629         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
630         """         """
631         return self.__shape         return self.__shape
# Line 544  class Symbol(object): Line 634  class Symbol(object):
634         """         """
635         the spatial dimension         the spatial dimension
636    
637         @return : the spatial dimension         @return: the spatial dimension
638         @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.
639         """         """
640         return self.__dim         return self.__dim
# Line 568  class Symbol(object): Line 658  class Symbol(object):
658         """         """
659         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.
660    
661         @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].
662         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
663         @rtype: C{list} of objects         @rtype: C{list} of objects
664         @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.
665         """         """
666         out=[]         out=[]
667         for a in self.getArgument():         for a in self.getArgument():
# Line 595  class Symbol(object): Line 685  class Symbol(object):
685            if isinstance(a,Symbol):            if isinstance(a,Symbol):
686               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
687            else:            else:
688                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
689                if len(s)>0:                if len(s)>0:
690                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
691                else:                else:
692                   out.append(a)                   out.append(a)
693         return out         return out
# Line 687  class Symbol(object): Line 777  class Symbol(object):
777         else:         else:
778            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
779            if len(s)>0:            if len(s)>0:
780               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
781            else:            else:
782               return 0.               return 0.
783    
# Line 695  class Symbol(object): Line 785  class Symbol(object):
785         """         """
786         returns -self.         returns -self.
787    
788         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
789         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
790         """         """
791         return self*(-1.)         return self*(-1.)
# Line 704  class Symbol(object): Line 794  class Symbol(object):
794         """         """
795         returns +self.         returns +self.
796    
797         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
798         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
799         """         """
800         return self*(1.)         return self*(1.)
801    
802     def __abs__(self):     def __abs__(self):
803         """         """
804         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
805         """         """
806         return Abs_Symbol(self)         return Abs_Symbol(self)
807    
# Line 721  class Symbol(object): Line 811  class Symbol(object):
811    
812         @param other: object to be added to this object         @param other: object to be added to this object
813         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
814         @return:  a S{Symbol} representing the sum of this object and C{other}         @return:  a L{Symbol} representing the sum of this object and C{other}
815         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
816         """         """
817         return add(self,other)         return add(self,other)
# Line 732  class Symbol(object): Line 822  class Symbol(object):
822    
823         @param other: object this object is added to         @param other: object this object is added to
824         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
825         @return: a S{Symbol} representing the sum of C{other} and this object object         @return: a L{Symbol} representing the sum of C{other} and this object object
826         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
827         """         """
828         return add(other,self)         return add(other,self)
# Line 743  class Symbol(object): Line 833  class Symbol(object):
833    
834         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
835         @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}.
836         @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
837         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
838         """         """
839         return add(self,-other)         return add(self,-other)
# Line 754  class Symbol(object): Line 844  class Symbol(object):
844    
845         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
846         @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}.
847         @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}.
848         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
849         """         """
850         return add(-self,other)         return add(-self,other)
# Line 765  class Symbol(object): Line 855  class Symbol(object):
855    
856         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
857         @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}.
858         @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}.
859         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
860         """         """
861         return mult(self,other)         return mult(self,other)
# Line 776  class Symbol(object): Line 866  class Symbol(object):
866    
867         @param other: object this object is multiplied with         @param other: object this object is multiplied with
868         @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}.
869         @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.
870         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
871         """         """
872         return mult(other,self)         return mult(other,self)
# Line 787  class Symbol(object): Line 877  class Symbol(object):
877    
878         @param other: object dividing this object         @param other: object dividing this object
879         @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}.
880         @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}
881         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
882         """         """
883         return quotient(self,other)         return quotient(self,other)
# Line 798  class Symbol(object): Line 888  class Symbol(object):
888    
889         @param other: object dividing this object         @param other: object dividing this object
890         @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}.
891         @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
892         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
893         """         """
894         return quotient(other,self)         return quotient(other,self)
# Line 809  class Symbol(object): Line 899  class Symbol(object):
899    
900         @param other: exponent         @param other: exponent
901         @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}.
902         @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}
903         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
904         """         """
905         return power(self,other)         return power(self,other)
# Line 820  class Symbol(object): Line 910  class Symbol(object):
910    
911         @param other: basis         @param other: basis
912         @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}.
913         @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
914         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
915         """         """
916         return power(other,self)         return power(other,self)
917    
918       def __getitem__(self,index):
919           """
920           returns the slice defined by index
921    
922           @param index: defines a
923           @type index: C{slice} or C{int} or a C{tuple} of them
924           @return: a L{Symbol} representing the slice defined by index
925           @rtype: L{DependendSymbol}
926           """
927           return GetSlice_Symbol(self,index)
928    
929  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
930     """     """
931     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.
932     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  
933        
934     Example:     Example::
935        
936     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
937     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
938     print u1==u2       print u1==u2
939     False       False
940        
941        but       but::
942    
943     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
944     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
945     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
946     print u1==u2, u1==u3       print u1==u2, u1==u3
947     True False       True False
948    
949     @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.
950     """     """
# Line 875  class DependendSymbol(Symbol): Line 976  class DependendSymbol(Symbol):
976  #=========================================================  #=========================================================
977  #  Unary operations prserving the shape  #  Unary operations prserving the shape
978  #========================================================  #========================================================
979    class GetSlice_Symbol(DependendSymbol):
980       """
981       L{Symbol} representing getting a slice for a L{Symbol}
982       """
983       def __init__(self,arg,index):
984          """
985          initialization of wherePositive L{Symbol} with argument arg
986          @param arg: argument
987          @type arg: L{Symbol}.
988          @param index: defines index
989          @type index: C{slice} or C{int} or a C{tuple} of them
990          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
991          @raises ValueError: if a step is given
992          """
993          if not isinstance(index,tuple): index=(index,)
994          if len(index)>arg.getRank():
995               raise IndexError,"GetSlice_Symbol: index out of range."
996          sh=()
997          index2=()
998          for i in range(len(index)):
999             ix=index[i]
1000             if isinstance(ix,int):
1001                if ix<0 or ix>=arg.getShape()[i]:
1002                   raise ValueError,"GetSlice_Symbol: index out of range."
1003                index2=index2+(ix,)
1004             else:
1005               if not ix.step==None:
1006                 raise ValueError,"GetSlice_Symbol: steping is not supported."
1007               if ix.start==None:
1008                  s=0
1009               else:
1010                  s=ix.start
1011               if ix.stop==None:
1012                  e=arg.getShape()[i]
1013               else:
1014                  e=ix.stop
1015                  if e>arg.getShape()[i]:
1016                     raise IndexError,"GetSlice_Symbol: index out of range."
1017               index2=index2+(slice(s,e),)
1018               if e>s:
1019                   sh=sh+(e-s,)
1020               elif s>e:
1021                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
1022          for i in range(len(index),arg.getRank()):
1023              index2=index2+(slice(0,arg.getShape()[i]),)
1024              sh=sh+(arg.getShape()[i],)
1025          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
1026    
1027       def getMyCode(self,argstrs,format="escript"):
1028          """
1029          returns a program code that can be used to evaluate the symbol.
1030    
1031          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1032          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1033          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1034          @type format: C{str}
1035          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1036          @rtype: C{str}
1037          @raise NotImplementedError: if the requested format is not available
1038          """
1039          if format=="escript" or format=="str"  or format=="text":
1040             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
1041          else:
1042             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
1043    
1044       def substitute(self,argvals):
1045          """
1046          assigns new values to symbols in the definition of the symbol.
1047          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1048    
1049          @param argvals: new values assigned to symbols
1050          @type argvals: C{dict} with keywords of type L{Symbol}.
1051          @return: result of the substitution process. Operations are executed as much as possible.
1052          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1053          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1054          """
1055          if argvals.has_key(self):
1056             arg=argvals[self]
1057             if self.isAppropriateValue(arg):
1058                return arg
1059             else:
1060                raise TypeError,"%s: new value is not appropriate."%str(self)
1061          else:
1062             args=self.getSubstitutedArguments(argvals)
1063             arg=args[0]
1064             index=args[1]
1065             return arg.__getitem__(index)
1066    
1067  def log10(arg):  def log10(arg):
1068     """     """
1069     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
1070    
1071     @param arg: argument     @param arg: argument
1072     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1073     @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.
1074     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1075     """     """
1076     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 903  def wherePositive(arg): Line 1092  def wherePositive(arg):
1092    
1093     @param arg: argument     @param arg: argument
1094     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1095     @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.
1096     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1097     """     """
1098     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1099        if arg.rank==0:        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1100           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1101             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))  
1102     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1103        return arg._wherePositive()        return arg._wherePositive()
1104     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 953  class WherePositive_Symbol(DependendSymb Line 1138  class WherePositive_Symbol(DependendSymb
1138        @type format: C{str}        @type format: C{str}
1139        @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.
1140        @rtype: C{str}        @rtype: C{str}
1141        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1142        """        """
1143        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1144            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 989  def whereNegative(arg): Line 1174  def whereNegative(arg):
1174    
1175     @param arg: argument     @param arg: argument
1176     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1177     @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.
1178     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1179     """     """
1180     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1181        if arg.rank==0:        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1182           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1183             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))  
1184     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1185        return arg._whereNegative()        return arg._whereNegative()
1186     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1039  class WhereNegative_Symbol(DependendSymb Line 1220  class WhereNegative_Symbol(DependendSymb
1220        @type format: C{str}        @type format: C{str}
1221        @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.
1222        @rtype: C{str}        @rtype: C{str}
1223        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1224        """        """
1225        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1226            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1075  def whereNonNegative(arg): Line 1256  def whereNonNegative(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.greater_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.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))  
1266     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1267        return arg._whereNonNegative()        return arg._whereNonNegative()
1268     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1109  def whereNonPositive(arg): Line 1286  def whereNonPositive(arg):
1286    
1287     @param arg: argument     @param arg: argument
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     @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.
1290     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1291     """     """
1292     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1293        if arg.rank==0:        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1294           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1295             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.  
1296     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1297        return arg._whereNonPositive()        return arg._whereNonPositive()
1298     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1145  def whereZero(arg,tol=0.): Line 1318  def whereZero(arg,tol=0.):
1318     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1319     @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.
1320     @type tol: C{float}     @type tol: C{float}
1321     @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.
1322     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1323     """     """
1324     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1325        if arg.rank==0:        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1326           if abs(arg)<=tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1327             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.  
1328     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1329        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1330     elif isinstance(arg,float):     elif isinstance(arg,float):
1331        if abs(arg)<=tol:        if abs(arg)<=tol:
1332          return 1.          return 1.
# Line 1198  class WhereZero_Symbol(DependendSymbol): Line 1364  class WhereZero_Symbol(DependendSymbol):
1364        @type format: C{str}        @type format: C{str}
1365        @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.
1366        @rtype: C{str}        @rtype: C{str}
1367        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1368        """        """
1369        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1370           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 1398  def whereNonZero(arg,tol=0.):
1398    
1399     @param arg: argument     @param arg: argument
1400     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @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.     @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.     @raises TypeError: if the type of the argument is not expected.
1403     """     """
1404     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1405        if arg.rank==0:        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1406          if abs(arg)>tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1407             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.  
1408     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1409        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1410     elif isinstance(arg,float):     elif isinstance(arg,float):
1411        if abs(arg)>tol:        if abs(arg)>tol:
1412          return 1.          return 1.
# Line 1263  def whereNonZero(arg,tol=0.): Line 1422  def whereNonZero(arg,tol=0.):
1422     else:     else:
1423        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1424    
1425    def erf(arg):
1426       """
1427       returns erf of argument arg
1428    
1429       @param arg: argument
1430       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1431       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1432       @raises TypeError: if the type of the argument is not expected.
1433       """
1434       if isinstance(arg,escript.Data):
1435          return arg._erf()
1436       else:
1437          raise TypeError,"erf: Unknown argument type."
1438    
1439  def sin(arg):  def sin(arg):
1440     """     """
1441     returns sine of argument arg     returns sine of argument arg
1442    
1443     @param arg: argument     @param arg: argument
1444     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1445     @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.
1446     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1447     """     """
1448     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1307  class Sin_Symbol(DependendSymbol): Line 1480  class Sin_Symbol(DependendSymbol):
1480        @type format: C{str}        @type format: C{str}
1481        @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.
1482        @rtype: C{str}        @rtype: C{str}
1483        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1484        """        """
1485        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1486            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1359  def cos(arg): Line 1532  def cos(arg):
1532    
1533     @param arg: argument     @param arg: argument
1534     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1535     @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.
1536     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1537     """     """
1538     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1397  class Cos_Symbol(DependendSymbol): Line 1570  class Cos_Symbol(DependendSymbol):
1570        @type format: C{str}        @type format: C{str}
1571        @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.
1572        @rtype: C{str}        @rtype: C{str}
1573        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1574        """        """
1575        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1576            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1449  def tan(arg): Line 1622  def tan(arg):
1622    
1623     @param arg: argument     @param arg: argument
1624     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1625     @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.
1626     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1627     """     """
1628     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1487  class Tan_Symbol(DependendSymbol): Line 1660  class Tan_Symbol(DependendSymbol):
1660        @type format: C{str}        @type format: C{str}
1661        @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.
1662        @rtype: C{str}        @rtype: C{str}
1663        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1664        """        """
1665        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1666            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1539  def asin(arg): Line 1712  def asin(arg):
1712    
1713     @param arg: argument     @param arg: argument
1714     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1715     @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.
1716     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1717     """     """
1718     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1577  class Asin_Symbol(DependendSymbol): Line 1750  class Asin_Symbol(DependendSymbol):
1750        @type format: C{str}        @type format: C{str}
1751        @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.
1752        @rtype: C{str}        @rtype: C{str}
1753        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1754        """        """
1755        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1756            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1629  def acos(arg): Line 1802  def acos(arg):
1802    
1803     @param arg: argument     @param arg: argument
1804     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1805     @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.
1806     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1807     """     """
1808     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1667  class Acos_Symbol(DependendSymbol): Line 1840  class Acos_Symbol(DependendSymbol):
1840        @type format: C{str}        @type format: C{str}
1841        @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.
1842        @rtype: C{str}        @rtype: C{str}
1843        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1844        """        """
1845        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1846            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1719  def atan(arg): Line 1892  def atan(arg):
1892    
1893     @param arg: argument     @param arg: argument
1894     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1895     @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.
1896     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1897     """     """
1898     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1757  class Atan_Symbol(DependendSymbol): Line 1930  class Atan_Symbol(DependendSymbol):
1930        @type format: C{str}        @type format: C{str}
1931        @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.
1932        @rtype: C{str}        @rtype: C{str}
1933        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1934        """        """
1935        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1936            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1809  def sinh(arg): Line 1982  def sinh(arg):
1982    
1983     @param arg: argument     @param arg: argument
1984     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1985     @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.
1986     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1987     """     """
1988     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1847  class Sinh_Symbol(DependendSymbol): Line 2020  class Sinh_Symbol(DependendSymbol):
2020        @type format: C{str}        @type format: C{str}
2021        @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.
2022        @rtype: C{str}        @rtype: C{str}
2023        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2024        """        """
2025        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2026            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1899  def cosh(arg): Line 2072  def cosh(arg):
2072    
2073     @param arg: argument     @param arg: argument
2074     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2075     @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.
2076     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2077     """     """
2078     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1937  class Cosh_Symbol(DependendSymbol): Line 2110  class Cosh_Symbol(DependendSymbol):
2110        @type format: C{str}        @type format: C{str}
2111        @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.
2112        @rtype: C{str}        @rtype: C{str}
2113        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2114        """        """
2115        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2116            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1989  def tanh(arg): Line 2162  def tanh(arg):
2162    
2163     @param arg: argument     @param arg: argument
2164     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2165     @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.
2166     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2167     """     """
2168     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2027  class Tanh_Symbol(DependendSymbol): Line 2200  class Tanh_Symbol(DependendSymbol):
2200        @type format: C{str}        @type format: C{str}
2201        @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.
2202        @rtype: C{str}        @rtype: C{str}
2203        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2204        """        """
2205        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2206            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2079  def asinh(arg): Line 2252  def asinh(arg):
2252    
2253     @param arg: argument     @param arg: argument
2254     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2255     @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.
2256     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2257     """     """
2258     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2117  class Asinh_Symbol(DependendSymbol): Line 2290  class Asinh_Symbol(DependendSymbol):
2290        @type format: C{str}        @type format: C{str}
2291        @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.
2292        @rtype: C{str}        @rtype: C{str}
2293        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2294        """        """
2295        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2296            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2169  def acosh(arg): Line 2342  def acosh(arg):
2342    
2343     @param arg: argument     @param arg: argument
2344     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2345     @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.
2346     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2347     """     """
2348     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2207  class Acosh_Symbol(DependendSymbol): Line 2380  class Acosh_Symbol(DependendSymbol):
2380        @type format: C{str}        @type format: C{str}
2381        @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.
2382        @rtype: C{str}        @rtype: C{str}
2383        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2384        """        """
2385        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2386            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2259  def atanh(arg): Line 2432  def atanh(arg):
2432    
2433     @param arg: argument     @param arg: argument
2434     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2435     @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.
2436     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2437     """     """
2438     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2297  class Atanh_Symbol(DependendSymbol): Line 2470  class Atanh_Symbol(DependendSymbol):
2470        @type format: C{str}        @type format: C{str}
2471        @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.
2472        @rtype: C{str}        @rtype: C{str}
2473        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2474        """        """
2475        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2476            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2349  def exp(arg): Line 2522  def exp(arg):
2522    
2523     @param arg: argument     @param arg: argument
2524     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2525     @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.
2526     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2527     """     """
2528     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2387  class Exp_Symbol(DependendSymbol): Line 2560  class Exp_Symbol(DependendSymbol):
2560        @type format: C{str}        @type format: C{str}
2561        @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.
2562        @rtype: C{str}        @rtype: C{str}
2563        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2564        """        """
2565        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2566            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2439  def sqrt(arg): Line 2612  def sqrt(arg):
2612    
2613     @param arg: argument     @param arg: argument
2614     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2615     @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.
2616     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2617     """     """
2618     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2477  class Sqrt_Symbol(DependendSymbol): Line 2650  class Sqrt_Symbol(DependendSymbol):
2650        @type format: C{str}        @type format: C{str}
2651        @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.
2652        @rtype: C{str}        @rtype: C{str}
2653        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2654        """        """
2655        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2656            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2529  def log(arg): Line 2702  def log(arg):
2702    
2703     @param arg: argument     @param arg: argument
2704     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2705     @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.
2706     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2707     """     """
2708     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2567  class Log_Symbol(DependendSymbol): Line 2740  class Log_Symbol(DependendSymbol):
2740        @type format: C{str}        @type format: C{str}
2741        @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.
2742        @rtype: C{str}        @rtype: C{str}
2743        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2744        """        """
2745        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2746            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2619  def sign(arg): Line 2792  def sign(arg):
2792    
2793     @param arg: argument     @param arg: argument
2794     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2795     @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.
2796     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2797     """     """
2798     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2667  class Abs_Symbol(DependendSymbol): Line 2840  class Abs_Symbol(DependendSymbol):
2840        @type format: C{str}        @type format: C{str}
2841        @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.
2842        @rtype: C{str}        @rtype: C{str}
2843        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2844        """        """
2845        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2846            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2719  def minval(arg): Line 2892  def minval(arg):
2892    
2893     @param arg: argument     @param arg: argument
2894     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2895     @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.
2896     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2897     """     """
2898     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2760  class Minval_Symbol(DependendSymbol): Line 2933  class Minval_Symbol(DependendSymbol):
2933        @type format: C{str}        @type format: C{str}
2934        @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.
2935        @rtype: C{str}        @rtype: C{str}
2936        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2937        """        """
2938        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2939            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2796  def maxval(arg): Line 2969  def maxval(arg):
2969    
2970     @param arg: argument     @param arg: argument
2971     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2972     @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.
2973     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2974     """     """
2975     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2837  class Maxval_Symbol(DependendSymbol): Line 3010  class Maxval_Symbol(DependendSymbol):
3010        @type format: C{str}        @type format: C{str}
3011        @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.
3012        @rtype: C{str}        @rtype: C{str}
3013        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3014        """        """
3015        if isinstance(argstrs,list):        if isinstance(argstrs,list):
3016            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2873  def length(arg): Line 3046  def length(arg):
3046    
3047     @param arg: argument     @param arg: argument
3048     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3049     @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.
3050     """     """
3051     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
3052    
3053    def trace(arg,axis_offset=0):
3054       """
3055       returns the trace of arg which the sum of arg[k,k] over k.
3056    
3057       @param arg: argument
3058       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3059       @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
3060                      C{axis_offset} and axis_offset+1 must be equal.
3061       @type axis_offset: C{int}
3062       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
3063       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3064       """
3065       if isinstance(arg,numarray.NumArray):
3066          sh=arg.shape
3067          if len(sh)<2:
3068            raise ValueError,"rank of argument must be greater than 1"
3069          if axis_offset<0 or axis_offset>len(sh)-2:
3070            raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
3071          s1=1
3072          for i in range(axis_offset): s1*=sh[i]
3073          s2=1
3074          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3075          if not sh[axis_offset] == sh[axis_offset+1]:
3076            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3077          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
3078          out=numarray.zeros([s1,s2],numarray.Float64)
3079          for i1 in range(s1):
3080            for i2 in range(s2):
3081                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
3082          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3083          return out
3084       elif isinstance(arg,escript.Data):
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          return arg._trace(axis_offset)
3093       elif isinstance(arg,float):
3094          raise TypeError,"illegal argument type float."
3095       elif isinstance(arg,int):
3096          raise TypeError,"illegal argument type int."
3097       elif isinstance(arg,Symbol):
3098          return Trace_Symbol(arg,axis_offset)
3099       else:
3100          raise TypeError,"Unknown argument type."
3101    
3102    class Trace_Symbol(DependendSymbol):
3103       """
3104       L{Symbol} representing the result of the trace function
3105       """
3106       def __init__(self,arg,axis_offset=0):
3107          """
3108          initialization of trace L{Symbol} with argument arg
3109          @param arg: argument of function
3110          @type arg: L{Symbol}.
3111          @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
3112                      C{axis_offset} and axis_offset+1 must be equal.
3113          @type axis_offset: C{int}
3114          """
3115          if arg.getRank()<2:
3116            raise ValueError,"rank of argument must be greater than 1"
3117          if axis_offset<0 or axis_offset>arg.getRank()-2:
3118            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3119          s=list(arg.getShape())        
3120          if not s[axis_offset] == s[axis_offset+1]:
3121            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3122          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3123    
3124       def getMyCode(self,argstrs,format="escript"):
3125          """
3126          returns a program code that can be used to evaluate the symbol.
3127    
3128          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3129          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3130          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3131          @type format: C{str}
3132          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3133          @rtype: C{str}
3134          @raise NotImplementedError: if the requested format is not available
3135          """
3136          if format=="escript" or format=="str"  or format=="text":
3137             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3138          else:
3139             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3140    
3141       def substitute(self,argvals):
3142          """
3143          assigns new values to symbols in the definition of the symbol.
3144          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3145    
3146          @param argvals: new values assigned to symbols
3147          @type argvals: C{dict} with keywords of type L{Symbol}.
3148          @return: result of the substitution process. Operations are executed as much as possible.
3149          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3150          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3151          """
3152          if argvals.has_key(self):
3153             arg=argvals[self]
3154             if self.isAppropriateValue(arg):
3155                return arg
3156             else:
3157                raise TypeError,"%s: new value is not appropriate."%str(self)
3158          else:
3159             arg=self.getSubstitutedArguments(argvals)
3160             return trace(arg[0],axis_offset=arg[1])
3161    
3162       def diff(self,arg):
3163          """
3164          differential of this object
3165    
3166          @param arg: the derivative is calculated with respect to arg
3167          @type arg: L{escript.Symbol}
3168          @return: derivative with respect to C{arg}
3169          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3170          """
3171          if arg==self:
3172             return identity(self.getShape())
3173          else:
3174             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3175    
3176    def transpose(arg,axis_offset=None):
3177       """
3178       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3179    
3180       @param arg: argument
3181       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3182       @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.
3183                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3184       @type axis_offset: C{int}
3185       @return: transpose of arg
3186       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3187       """
3188       if isinstance(arg,numarray.NumArray):
3189          if axis_offset==None: axis_offset=int(arg.rank/2)
3190          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3191       elif isinstance(arg,escript.Data):
3192          r=arg.getRank()
3193          if axis_offset==None: axis_offset=int(r/2)
3194          if axis_offset<0 or axis_offset>r:
3195            raise ValueError,"axis_offset must be between 0 and %s"%r
3196          return arg._transpose(axis_offset)
3197       elif isinstance(arg,float):
3198          if not ( axis_offset==0 or axis_offset==None):
3199            raise ValueError,"axis_offset must be 0 for float argument"
3200          return arg
3201       elif isinstance(arg,int):
3202          if not ( axis_offset==0 or axis_offset==None):
3203            raise ValueError,"axis_offset must be 0 for int argument"
3204          return float(arg)
3205       elif isinstance(arg,Symbol):
3206          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3207          return Transpose_Symbol(arg,axis_offset)
3208       else:
3209          raise TypeError,"Unknown argument type."
3210    
3211    class Transpose_Symbol(DependendSymbol):
3212       """
3213       L{Symbol} representing the result of the transpose function
3214       """
3215       def __init__(self,arg,axis_offset=None):
3216          """
3217          initialization of transpose L{Symbol} with argument arg
3218    
3219          @param arg: argument of function
3220          @type arg: L{Symbol}.
3221          @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.
3222                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3223          @type axis_offset: C{int}
3224          """
3225          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3226          if axis_offset<0 or axis_offset>arg.getRank():
3227            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3228          s=arg.getShape()
3229          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3230    
3231       def getMyCode(self,argstrs,format="escript"):
3232          """
3233          returns a program code that can be used to evaluate the symbol.
3234    
3235          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3236          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3237          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3238          @type format: C{str}
3239          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3240          @rtype: C{str}
3241          @raise NotImplementedError: if the requested format is not available
3242          """
3243          if format=="escript" or format=="str"  or format=="text":
3244             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3245          else:
3246             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3247    
3248       def substitute(self,argvals):
3249          """
3250          assigns new values to symbols in the definition of the symbol.
3251          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3252    
3253          @param argvals: new values assigned to symbols
3254          @type argvals: C{dict} with keywords of type L{Symbol}.
3255          @return: result of the substitution process. Operations are executed as much as possible.
3256          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3257          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3258          """
3259          if argvals.has_key(self):
3260             arg=argvals[self]
3261             if self.isAppropriateValue(arg):
3262                return arg
3263             else:
3264                raise TypeError,"%s: new value is not appropriate."%str(self)
3265          else:
3266             arg=self.getSubstitutedArguments(argvals)
3267             return transpose(arg[0],axis_offset=arg[1])
3268    
3269       def diff(self,arg):
3270          """
3271          differential of this object
3272    
3273          @param arg: the derivative is calculated with respect to arg
3274          @type arg: L{escript.Symbol}
3275          @return: derivative with respect to C{arg}
3276          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3277          """
3278          if arg==self:
3279             return identity(self.getShape())
3280          else:
3281             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3282    
3283    def swap_axes(arg,axis0=0,axis1=1):
3284       """
3285       returns the swap of arg by swaping the components axis0 and axis1
3286    
3287       @param arg: argument
3288       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3289       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3290       @type axis0: C{int}
3291       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3292       @type axis1: C{int}
3293       @return: C{arg} with swaped components
3294       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3295       """
3296       if axis0 > axis1:
3297          axis0,axis1=axis1,axis0
3298       if isinstance(arg,numarray.NumArray):
3299          return numarray.swapaxes(arg,axis0,axis1)
3300       elif isinstance(arg,escript.Data):
3301          return arg._swap_axes(axis0,axis1)
3302       elif isinstance(arg,float):
3303          raise TyepError,"float argument is not supported."
3304       elif isinstance(arg,int):
3305          raise TyepError,"int argument is not supported."
3306       elif isinstance(arg,Symbol):
3307          return SwapAxes_Symbol(arg,axis0,axis1)
3308       else:
3309          raise TypeError,"Unknown argument type."
3310    
3311    class SwapAxes_Symbol(DependendSymbol):
3312       """
3313       L{Symbol} representing the result of the swap function
3314       """
3315       def __init__(self,arg,axis0=0,axis1=1):
3316          """
3317          initialization of swap L{Symbol} with argument arg
3318    
3319          @param arg: argument
3320          @type arg: L{Symbol}.
3321          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3322          @type axis0: C{int}
3323          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3324          @type axis1: C{int}
3325          """
3326          if arg.getRank()<2:
3327             raise ValueError,"argument must have at least rank 2."
3328          if axis0<0 or axis0>arg.getRank()-1:
3329             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3330          if axis1<0 or axis1>arg.getRank()-1:
3331             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3332          if axis0 == axis1:
3333             raise ValueError,"axis indices must be different."
3334          if axis0 > axis1:
3335             axis0,axis1=axis1,axis0
3336          s=arg.getShape()
3337          s_out=[]
3338          for i in range(len(s)):
3339             if i == axis0:
3340                s_out.append(s[axis1])
3341             elif i == axis1:
3342                s_out.append(s[axis0])
3343             else:
3344                s_out.append(s[i])
3345          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3346    
3347       def getMyCode(self,argstrs,format="escript"):
3348          """
3349          returns a program code that can be used to evaluate the symbol.
3350    
3351          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3352          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3353          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3354          @type format: C{str}
3355          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3356          @rtype: C{str}
3357          @raise NotImplementedError: if the requested format is not available
3358          """
3359          if format=="escript" or format=="str"  or format=="text":
3360             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3361          else:
3362             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3363    
3364       def substitute(self,argvals):
3365          """
3366          assigns new values to symbols in the definition of the symbol.
3367          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3368    
3369          @param argvals: new values assigned to symbols
3370          @type argvals: C{dict} with keywords of type L{Symbol}.
3371          @return: result of the substitution process. Operations are executed as much as possible.
3372          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3373          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3374          """
3375          if argvals.has_key(self):
3376             arg=argvals[self]
3377             if self.isAppropriateValue(arg):
3378                return arg
3379             else:
3380                raise TypeError,"%s: new value is not appropriate."%str(self)
3381          else:
3382             arg=self.getSubstitutedArguments(argvals)
3383             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3384    
3385       def diff(self,arg):
3386          """
3387          differential of this object
3388    
3389          @param arg: the derivative is calculated with respect to arg
3390          @type arg: L{escript.Symbol}
3391          @return: derivative with respect to C{arg}
3392          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3393          """
3394          if arg==self:
3395             return identity(self.getShape())
3396          else:
3397             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3398    
3399    def symmetric(arg):
3400        """
3401        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3402    
3403        @param arg: square matrix. Must have rank 2 or 4 and be square.
3404        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3405        @return: symmetric part of arg
3406        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3407        """
3408        if isinstance(arg,numarray.NumArray):
3409          if arg.rank==2:
3410            if not (arg.shape[0]==arg.shape[1]):
3411               raise ValueError,"argument must be square."
3412          elif arg.rank==4:
3413            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3414               raise ValueError,"argument must be square."
3415          else:
3416            raise ValueError,"rank 2 or 4 is required."
3417          return (arg+transpose(arg))/2
3418        elif isinstance(arg,escript.Data):
3419          if arg.getRank()==2:
3420            if not (arg.getShape()[0]==arg.getShape()[1]):
3421               raise ValueError,"argument must be square."
3422            return arg._symmetric()
3423          elif arg.getRank()==4:
3424            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3425               raise ValueError,"argument must be square."
3426            return arg._symmetric()
3427          else:
3428            raise ValueError,"rank 2 or 4 is required."
3429        elif isinstance(arg,float):
3430          return arg
3431        elif isinstance(arg,int):
3432          return float(arg)
3433        elif isinstance(arg,Symbol):
3434          if arg.getRank()==2:
3435            if not (arg.getShape()[0]==arg.getShape()[1]):
3436               raise ValueError,"argument must be square."
3437          elif arg.getRank()==4:
3438            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3439               raise ValueError,"argument must be square."
3440          else:
3441            raise ValueError,"rank 2 or 4 is required."
3442          return (arg+transpose(arg))/2
3443        else:
3444          raise TypeError,"symmetric: Unknown argument type."
3445    
3446    def nonsymmetric(arg):
3447        """
3448        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3449    
3450        @param arg: square matrix. Must have rank 2 or 4 and be square.
3451        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3452        @return: nonsymmetric part of arg
3453        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3454        """
3455        if isinstance(arg,numarray.NumArray):
3456          if arg.rank==2:
3457            if not (arg.shape[0]==arg.shape[1]):
3458               raise ValueError,"nonsymmetric: argument must be square."
3459          elif arg.rank==4:
3460            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3461               raise ValueError,"nonsymmetric: argument must be square."
3462          else:
3463            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3464          return (arg-transpose(arg))/2
3465        elif isinstance(arg,escript.Data):
3466          if arg.getRank()==2:
3467            if not (arg.getShape()[0]==arg.getShape()[1]):
3468               raise ValueError,"argument must be square."
3469            return arg._nonsymmetric()
3470          elif arg.getRank()==4:
3471            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3472               raise ValueError,"argument must be square."
3473            return arg._nonsymmetric()
3474          else:
3475            raise ValueError,"rank 2 or 4 is required."
3476        elif isinstance(arg,float):
3477          return arg
3478        elif isinstance(arg,int):
3479          return float(arg)
3480        elif isinstance(arg,Symbol):
3481          if arg.getRank()==2:
3482            if not (arg.getShape()[0]==arg.getShape()[1]):
3483               raise ValueError,"nonsymmetric: argument must be square."
3484          elif arg.getRank()==4:
3485            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3486               raise ValueError,"nonsymmetric: argument must be square."
3487          else:
3488            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3489          return (arg-transpose(arg))/2
3490        else:
3491          raise TypeError,"nonsymmetric: Unknown argument type."
3492    
3493    def inverse(arg):
3494        """
3495        returns the inverse of the square matrix arg.
3496    
3497        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3498        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3499        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3500        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3501        @note: for L{escript.Data} objects the dimension is restricted to 3.
3502        """
3503        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3504        if isinstance(arg,numarray.NumArray):
3505          return numarray.linear_algebra.inverse(arg)
3506        elif isinstance(arg,escript.Data):
3507          return escript_inverse(arg)
3508        elif isinstance(arg,float):
3509          return 1./arg
3510        elif isinstance(arg,int):
3511          return 1./float(arg)
3512        elif isinstance(arg,Symbol):
3513          return Inverse_Symbol(arg)
3514        else:
3515          raise TypeError,"inverse: Unknown argument type."
3516    
3517    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3518          "arg is a Data objects!!!"
3519          if not arg.getRank()==2:
3520            raise ValueError,"escript_inverse: argument must have rank 2"
3521          s=arg.getShape()      
3522          if not s[0] == s[1]:
3523            raise ValueError,"escript_inverse: argument must be a square matrix."
3524          out=escript.Data(0.,s,arg.getFunctionSpace())
3525          if s[0]==1:
3526              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3527                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3528              out[0,0]=1./arg[0,0]
3529          elif s[0]==2:
3530              A11=arg[0,0]
3531              A12=arg[0,1]
3532              A21=arg[1,0]
3533              A22=arg[1,1]
3534              D = A11*A22-A12*A21
3535              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3536                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3537              D=1./D
3538              out[0,0]= A22*D
3539              out[1,0]=-A21*D
3540              out[0,1]=-A12*D
3541              out[1,1]= A11*D
3542          elif s[0]==3:
3543              A11=arg[0,0]
3544              A21=arg[1,0]
3545              A31=arg[2,0]
3546              A12=arg[0,1]
3547              A22=arg[1,1]
3548              A32=arg[2,1]
3549              A13=arg[0,2]
3550              A23=arg[1,2]
3551              A33=arg[2,2]
3552              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3553              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3554                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3555              D=1./D
3556              out[0,0]=(A22*A33-A23*A32)*D
3557              out[1,0]=(A31*A23-A21*A33)*D
3558              out[2,0]=(A21*A32-A31*A22)*D
3559              out[0,1]=(A13*A32-A12*A33)*D
3560              out[1,1]=(A11*A33-A31*A13)*D
3561              out[2,1]=(A12*A31-A11*A32)*D
3562              out[0,2]=(A12*A23-A13*A22)*D
3563              out[1,2]=(A13*A21-A11*A23)*D
3564              out[2,2]=(A11*A22-A12*A21)*D
3565          else:
3566             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3567          return out
3568    
3569    class Inverse_Symbol(DependendSymbol):
3570       """
3571       L{Symbol} representing the result of the inverse function
3572       """
3573       def __init__(self,arg):
3574          """
3575          initialization of inverse L{Symbol} with argument arg
3576          @param arg: argument of function
3577          @type arg: L{Symbol}.
3578          """
3579          if not arg.getRank()==2:
3580            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3581          s=arg.getShape()
3582          if not s[0] == s[1]:
3583            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3584          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3585    
3586       def getMyCode(self,argstrs,format="escript"):
3587          """
3588          returns a program code that can be used to evaluate the symbol.
3589    
3590          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3591          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3592          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3593          @type format: C{str}
3594          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3595          @rtype: C{str}
3596          @raise NotImplementedError: if the requested format is not available
3597          """
3598          if format=="escript" or format=="str"  or format=="text":
3599             return "inverse(%s)"%argstrs[0]
3600          else:
3601             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3602    
3603       def substitute(self,argvals):
3604          """
3605          assigns new values to symbols in the definition of the symbol.
3606          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3607    
3608          @param argvals: new values assigned to symbols
3609          @type argvals: C{dict} with keywords of type L{Symbol}.
3610          @return: result of the substitution process. Operations are executed as much as possible.
3611          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3612          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3613          """
3614          if argvals.has_key(self):
3615             arg=argvals[self]
3616             if self.isAppropriateValue(arg):
3617                return arg
3618             else:
3619                raise TypeError,"%s: new value is not appropriate."%str(self)
3620          else:
3621             arg=self.getSubstitutedArguments(argvals)
3622             return inverse(arg[0])
3623    
3624       def diff(self,arg):
3625          """
3626          differential of this object
3627    
3628          @param arg: the derivative is calculated with respect to arg
3629          @type arg: L{escript.Symbol}
3630          @return: derivative with respect to C{arg}
3631          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3632          """
3633          if arg==self:
3634             return identity(self.getShape())
3635          else:
3636             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3637    
3638    def eigenvalues(arg):
3639        """
3640        returns the eigenvalues of the square matrix arg.
3641    
3642        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3643                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3644        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3645        @return: the eigenvalues in increasing order.
3646        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3647        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3648        """
3649        if isinstance(arg,numarray.NumArray):
3650          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3651          out.sort()
3652          return out
3653        elif isinstance(arg,escript.Data):
3654          return arg._eigenvalues()
3655        elif isinstance(arg,Symbol):
3656          if not arg.getRank()==2:
3657            raise ValueError,"eigenvalues: argument must have rank 2"
3658          s=arg.getShape()      
3659          if not s[0] == s[1]:
3660            raise ValueError,"eigenvalues: argument must be a square matrix."
3661          if s[0]==1:
3662              return arg[0]
3663          elif s[0]==2:
3664              arg1=symmetric(arg)
3665              A11=arg1[0,0]
3666              A12=arg1[0,1]
3667              A22=arg1[1,1]
3668              trA=(A11+A22)/2.
3669              A11-=trA
3670              A22-=trA
3671              s=sqrt(A12**2-A11*A22)
3672              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3673          elif s[0]==3:
3674              arg1=symmetric(arg)
3675              A11=arg1[0,0]
3676              A12=arg1[0,1]
3677              A22=arg1[1,1]
3678              A13=arg1[0,2]
3679              A23=arg1[1,2]
3680              A33=arg1[2,2]
3681              trA=(A11+A22+A33)/3.
3682              A11-=trA
3683              A22-=trA
3684              A33-=trA
3685              A13_2=A13**2
3686              A23_2=A23**2
3687              A12_2=A12**2
3688              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3689              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3690              sq_p=sqrt(p/3.)
3691              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
3692              sq_p*=2.
3693              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3694               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3695               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3696              return trA+sq_p*f
3697          else:
3698             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3699        elif isinstance(arg,float):
3700          return arg
3701        elif isinstance(arg,int):
3702          return float(arg)
3703        else:
3704          raise TypeError,"eigenvalues: Unknown argument type."
3705    
3706    def eigenvalues_and_eigenvectors(arg):
3707        """
3708        returns the eigenvalues and eigenvectors of the square matrix arg.
3709    
3710        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3711                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3712        @type arg: L{escript.Data}
3713        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3714                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3715                 the eigenvector coresponding to the i-th eigenvalue.
3716        @rtype: L{tuple} of L{escript.Data}.
3717        @note: The dimension is restricted to 3.
3718        """
3719        if isinstance(arg,numarray.NumArray):
3720          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3721        elif isinstance(arg,escript.Data):
3722          return arg._eigenvalues_and_eigenvectors()
3723        elif isinstance(arg,Symbol):
3724          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3725        elif isinstance(arg,float):
3726          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3727        elif isinstance(arg,int):
3728          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3729        else:
3730          raise TypeError,"eigenvalues: Unknown argument type."
3731  #=======================================================  #=======================================================
3732  #  Binary operations:  #  Binary operations:
3733  #=======================================================  #=======================================================
# Line 2920  class Add_Symbol(DependendSymbol): Line 3771  class Add_Symbol(DependendSymbol):
3771         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3772         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3773         """         """
3774         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3775         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3776         if not sh0==sh1:         if not sh0==sh1:
3777            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3778         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 3787  class Add_Symbol(DependendSymbol):
3787        @type format: C{str}        @type format: C{str}
3788        @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.
3789        @rtype: C{str}        @rtype: C{str}
3790        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3791        """        """
3792        if format=="str" or format=="text":        if format=="str" or format=="text":
3793           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 2964  class Add_Symbol(DependendSymbol): Line 3815  class Add_Symbol(DependendSymbol):
3815              raise TypeError,"%s: new value is not appropriate."%str(self)              raise TypeError,"%s: new value is not appropriate."%str(self)
3816        else:        else:
3817           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
3818           return add(args[0],args[1])           out=add(args[0],args[1])
3819             return out
3820    
3821     def diff(self,arg):     def diff(self,arg):
3822        """        """
# Line 2995  def mult(arg0,arg1): Line 3847  def mult(arg0,arg1):
3847         """         """
3848         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3849         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3850            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3851         else:         else:
3852            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3853                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3019  class Mult_Symbol(DependendSymbol): Line 3871  class Mult_Symbol(DependendSymbol):
3871         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3872         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3873         """         """
3874         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3875         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3876         if not sh0==sh1:         if not sh0==sh1:
3877            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3878         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 3887  class Mult_Symbol(DependendSymbol):
3887        @type format: C{str}        @type format: C{str}
3888        @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.
3889        @rtype: C{str}        @rtype: C{str}
3890        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3891        """        """
3892        if format=="str" or format=="text":        if format=="str" or format=="text":
3893           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3095  def quotient(arg0,arg1): Line 3947  def quotient(arg0,arg1):
3947         """         """
3948         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3949         if testForZero(args[0]):         if testForZero(args[0]):
3950            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3951         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3952            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3953               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3124  class Quotient_Symbol(DependendSymbol): Line 3976  class Quotient_Symbol(DependendSymbol):
3976         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3977         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3978         """         """
3979         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3980         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3981         if not sh0==sh1:         if not sh0==sh1:
3982            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3983         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 3992  class Quotient_Symbol(DependendSymbol):
3992        @type format: C{str}        @type format: C{str}
3993        @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.
3994        @rtype: C{str}        @rtype: C{str}
3995        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3996        """        """
3997        if format=="str" or format=="text":        if format=="str" or format=="text":
3998           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3201  def power(arg0,arg1): Line 4053  def power(arg0,arg1):
4053         """         """
4054         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4055         if testForZero(args[0]):         if testForZero(args[0]):
4056            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4057         elif testForZero(args[1]):         elif testForZero(args[1]):
4058            return numarray.ones(args[0],numarray.Float)            return numarray.ones(getShape(args[1]),numarray.Float64)
4059         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4060            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4061         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 4078  class Power_Symbol(DependendSymbol):
4078         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4079         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4080         """         """
4081         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4082         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4083         if not sh0==sh1:         if not sh0==sh1:
4084            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4085         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3244  class Power_Symbol(DependendSymbol): Line 4096  class Power_Symbol(DependendSymbol):
4096        @type format: C{str}        @type format: C{str}
4097        @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.
4098        @rtype: C{str}        @rtype: C{str}
4099        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4100        """        """
4101        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4102           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3304  def maximum(*args): Line 4156  def maximum(*args):
4156         if out==None:         if out==None:
4157            out=a            out=a
4158         else:         else:
4159            m=whereNegative(out-a)            diff=add(a,-out)
4160            out=m*a+(1.-m)*out            out=add(out,mult(wherePositive(diff),diff))
4161      return out      return out
4162        
4163  def minimum(*arg):  def minimum(*args):
4164      """      """
4165      the minimum over arguments args      the minimum over arguments args
4166    
# Line 3322  def minimum(*arg): Line 4174  def minimum(*arg):
4174         if out==None:         if out==None:
4175            out=a            out=a
4176         else:         else:
4177            m=whereNegative(out-a)            diff=add(a,-out)
4178            out=m*out+(1.-m)*a            out=add(out,mult(whereNegative(diff),diff))
4179      return out      return out
4180    
4181    def clip(arg,minval=None,maxval=None):
4182        """
4183        cuts the values of arg between minval and maxval
4184    
4185        @param arg: argument
4186        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4187        @param minval: lower range. If None no lower range is applied
4188        @type minval: C{float} or C{None}
4189        @param maxval: upper range. If None no upper range is applied
4190        @type maxval: C{float} or C{None}
4191        @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4192                 less then maxval are unchanged.
4193        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4194        @raise ValueError: if minval>maxval
4195        """
4196        if not minval==None and not maxval==None:
4197           if minval>maxval:
4198              raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4199        if minval == None:
4200            tmp=arg
4201        else:
4202            tmp=maximum(minval,arg)
4203        if maxval == None:
4204            return tmp
4205        else:
4206            return minimum(tmp,maxval)
4207    
4208        
4209  def inner(arg0,arg1):  def inner(arg0,arg1):
4210      """      """
# Line 3340  def inner(arg0,arg1): Line 4220  def inner(arg0,arg1):
4220      @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}
4221      @param arg1: second argument      @param arg1: second argument
4222      @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}
4223      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4224      @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
4225      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4226      """      """
4227      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4228      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4229      if not sh0==sh1:      if not sh0==sh1:
4230          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4231      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4232    
4233    def outer(arg0,arg1):
4234        """
4235        the outer product of the two argument:
4236    
4237        out[t,s]=arg0[t]*arg1[s]
4238    
4239        where
4240    
4241            - s runs through arg0.Shape
4242            - t runs through arg1.Shape
4243    
4244        @param arg0: first argument
4245        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4246        @param arg1: second argument
4247        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4248        @return: the outer product of arg0 and arg1 at each data point
4249        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4250        """
4251        return generalTensorProduct(arg0,arg1,axis_offset=0)
4252    
4253  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4254      """      """
4255        see L{matrix_mult}
4256        """
4257        return matrix_mult(arg0,arg1)
4258    
4259    def matrix_mult(arg0,arg1):
4260        """
4261      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4262    
4263      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4264    
4265            or      or
4266    
4267      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4268    
4269      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.
4270    
4271      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4272      @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 4276  def matrixmult(arg0,arg1):
4276      @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
4277      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4278      """      """
4279      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4280      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4281      if not len(sh0)==2 :      if not len(sh0)==2 :
4282          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4283      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4284          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4285      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4286    
4287  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4288      """      """
4289      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  
4290      """      """
4291      return generalTensorProduct(arg0,arg1,offset=0)      return tensor_mult(arg0,arg1)
4292    
4293    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4294      """      """
4295      the tensor product of the two argument:      the tensor product of the two argument:
   
4296            
4297      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4298    
4299      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4300    
4301                   or      or
4302    
4303      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4304    
# Line 3415  def tensormult(arg0,arg1): Line 4307  def tensormult(arg0,arg1):
4307    
4308      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]
4309                                
4310                   or      or
4311    
4312      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]
4313    
4314                   or      or
4315    
4316      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]
4317    
4318      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  
4319      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.
4320    
4321      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4322      @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 4325  def tensormult(arg0,arg1):
4325      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4326      @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
4327      """      """
4328      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4329      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4330      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4331         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4332      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):
4333         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4334      else:      else:
4335          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4336    
4337  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4338      """      """
4339      generalized tensor product      generalized tensor product
4340    
4341      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4342    
4343      where s runs through arg0.Shape[:arg0.Rank-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[offset:]  
4344    
4345      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]
4346      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4347            - t runs through arg1.Shape[axis_offset:]
4348    
4349      @param arg0: first argument      @param arg0: first argument
4350      @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}
4351      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4352      @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}
4353      @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.
4354      @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 4358  def generalTensorProduct(arg0,arg1,offse
4358      # 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
4359      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4360         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4361             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4362         else:         else:
4363             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4364                 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)
4365             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4366             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4367             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4368             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4369             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4370             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4371             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4372             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4373             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4374             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4375             for i0 in range(d0):             for i0 in range(d0):
4376                      for i1 in range(d1):                      for i1 in range(d1):
4377                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4378             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4379             return out             return out
4380      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4381         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4382             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4383         else:         else:
4384             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)
4385      else:            else:      
4386         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4387                                    
4388  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4389     """     """
4390     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4391     """     """
4392     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4393         """         """
4394         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4395    
4396         @param arg0: numerator         @param arg0: first argument
4397         @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}.
4398         @param arg1: denominator         @param arg1: second argument
4399         @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}.
4400         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4401         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4402         """         """
4403         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4404         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4405         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4406         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4407         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4408         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4409         if not sh01==sh10:         if not sh01==sh10:
4410             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)
4411         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])
4412    
4413     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4414        """        """
# Line 3529  class GeneralTensorProduct_Symbol(Depend Line 4420  class GeneralTensorProduct_Symbol(Depend
4420        @type format: C{str}        @type format: C{str}
4421        @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.
4422        @rtype: C{str}        @rtype: C{str}
4423        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4424        """        """
4425        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4426           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])
4427        else:        else:
4428           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)
4429    
# Line 3557  class GeneralTensorProduct_Symbol(Depend Line 4448  class GeneralTensorProduct_Symbol(Depend
4448           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4449           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4450    
4451  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4452      "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!!!"
4453      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4454      shape0=arg0.getShape()[:arg0.getRank()-offset]  
4455      shape01=arg0.getShape()[arg0.getRank()-offset:]  def transposed_matrix_mult(arg0,arg1):
4456      shape10=arg1.getShape()[:offset]      """
4457      shape1=arg1.getShape()[offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4458      if not shape01==shape10:  
4459          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]
4460    
4461      # create return value:      or
4462      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())  
4463      #      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4464      s0=[[]]  
4465      for k in shape0:      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4466            s=[]  
4467            for j in s0:      The first dimension of arg0 and arg1 must match.
4468                  for i in range(k): s.append(j+[slice(i,i)])  
4469            s0=s      @param arg0: first argument of rank 2
4470      s1=[[]]      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4471      for k in shape1:      @param arg1: second argument of at least rank 1
4472            s=[]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4473            for j in s1:      @return: the product of the transposed of arg0 and arg1 at each data point
4474                  for i in range(k): s.append(j+[slice(i,i)])      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4475            s1=s      @raise ValueError: if the shapes of the arguments are not appropriate
4476      s01=[[]]      """
4477      for k in shape01:      sh0=getShape(arg0)
4478            s=[]      sh1=getShape(arg1)
4479            for j in s01:      if not len(sh0)==2 :
4480                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"first argument must have rank 2"
4481            s01=s      if not len(sh1)==2 and not len(sh1)==1:
4482            raise ValueError,"second argument must have rank 1 or 2"
4483      for i0 in s0:      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4484         for i1 in s1:  
4485           s=escript.Scalar(0.,arg0.getFunctionSpace())  def transposed_tensor_mult(arg0,arg1):
4486           for i01 in s01:      """
4487              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      the tensor product of the transposed of the first and the second argument
4488           out.__setitem__(tuple(i0+i1),s)      
4489      return out      for arg0 of rank 2 this is
4490    
4491        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4492    
4493        or
4494    
4495        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4496    
4497      
4498        and for arg0 of rank 4 this is
4499    
4500        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4501                  
4502        or
4503    
4504        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4505    
4506        or
4507    
4508        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4509    
4510        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4511        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4512    
4513        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4514    
4515        @param arg0: first argument of rank 2 or 4
4516        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4517        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4518        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4519        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4520        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4521        """
4522        sh0=getShape(arg0)
4523        sh1=getShape(arg1)
4524        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4525           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4526        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4527           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4528        else:
4529            raise ValueError,"first argument must have rank 2 or 4"
4530    
4531    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4532        """
4533        generalized tensor product of transposed of arg0 and arg1:
4534    
4535        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4536    
4537        where
4538    
4539            - s runs through arg0.Shape[axis_offset:]
4540            - r runs trough arg0.Shape[:axis_offset]
4541            - t runs through arg1.Shape[axis_offset:]
4542    
4543        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4544        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4545    
4546        @param arg0: first argument
4547        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4548        @param arg1: second argument
4549        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4550        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4551        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4552        """
4553        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4554        arg0,arg1=matchType(arg0,arg1)
4555        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4556        if isinstance(arg0,numarray.NumArray):
4557           if isinstance(arg1,Symbol):
4558               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4559           else:
4560               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4561                   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)
4562               arg0_c=arg0.copy()
4563               arg1_c=arg1.copy()
4564               sh0,sh1=arg0.shape,arg1.shape
4565               d0,d1,d01=1,1,1
4566               for i in sh0[axis_offset:]: d0*=i
4567               for i in sh1[axis_offset:]: d1*=i
4568               for i in sh0[:axis_offset]: d01*=i
4569               arg0_c.resize((d01,d0))
4570               arg1_c.resize((d01,d1))
4571               out=numarray.zeros((d0,d1),numarray.Float64)
4572               for i0 in range(d0):
4573                        for i1 in range(d1):
4574                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4575               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4576               return out
4577        elif isinstance(arg0,escript.Data):
4578           if isinstance(arg1,Symbol):
4579               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4580           else:
4581               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4582        else:      
4583           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4584                    
4585    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4586       """
4587       Symbol representing the general tensor product of the transposed of arg0 and arg1
4588       """
4589       def __init__(self,arg0,arg1,axis_offset=0):
4590           """
4591           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4592    
4593           @param arg0: first argument
4594           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4595           @param arg1: second argument
4596           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4597           @raise ValueError: inconsistent dimensions of arguments.
4598           @note: if both arguments have a spatial dimension, they must equal.
4599           """
4600           sh_arg0=getShape(arg0)
4601           sh_arg1=getShape(arg1)
4602           sh01=sh_arg0[:axis_offset]
4603           sh10=sh_arg1[:axis_offset]
4604           sh0=sh_arg0[axis_offset:]
4605           sh1=sh_arg1[axis_offset:]
4606           if not sh01==sh10:
4607               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)
4608           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4609    
4610       def getMyCode(self,argstrs,format="escript"):
4611          """
4612          returns a program code that can be used to evaluate the symbol.
4613    
4614          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4615          @type argstrs: C{list} of length 2 of C{str}.
4616          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4617          @type format: C{str}
4618          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4619          @rtype: C{str}
4620          @raise NotImplementedError: if the requested format is not available
4621          """
4622          if format=="escript" or format=="str" or format=="text":
4623             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4624          else:
4625             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4626    
4627       def substitute(self,argvals):
4628          """
4629          assigns new values to symbols in the definition of the symbol.
4630          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4631    
4632          @param argvals: new values assigned to symbols
4633          @type argvals: C{dict} with keywords of type L{Symbol}.
4634          @return: result of the substitution process. Operations are executed as much as possible.
4635          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4636          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4637          """
4638          if argvals.has_key(self):
4639             arg=argvals[self]
4640             if self.isAppropriateValue(arg):
4641                return arg
4642             else:
4643                raise TypeError,"%s: new value is not appropriate."%str(self)
4644          else:
4645             args=self.getSubstitutedArguments(argvals)
4646             return generalTransposedTensorProduct(args[0],args[1],args[2])
4647    
4648    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4649        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4650        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4651    
4652    def matrix_transposed_mult(arg0,arg1):
4653        """
4654        matrix-transposed(matrix) product of the two argument:
4655    
4656        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4657    
4658        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4659    
4660        The last dimensions of arg0 and arg1 must match.
4661    
4662        @param arg0: first argument of rank 2
4663        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4664        @param arg1: second argument of rank 2
4665        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4666        @return: the product of arg0 and the transposed of arg1 at each data point
4667        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4668        @raise ValueError: if the shapes of the arguments are not appropriate
4669        """
4670        sh0=getShape(arg0)
4671        sh1=getShape(arg1)
4672        if not len(sh0)==2 :
4673            raise ValueError,"first argument must have rank 2"
4674        if not len(sh1)==2 and not len(sh1)==1:
4675            raise ValueError,"second argument must have rank 1 or 2"
4676        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4677    
4678    def tensor_transposed_mult(arg0,arg1):
4679        """
4680        the tensor product of the first and the transpose of the second argument
4681        
4682        for arg0 of rank 2 this is
4683    
4684        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4685    
4686        and for arg0 of rank 4 this is
4687    
4688        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4689                  
4690        or
4691    
4692        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4693    
4694        In the first case the the second dimension of arg0 and arg1 must match and  
4695        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4696    
4697        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4698    
4699        @param arg0: first argument of rank 2 or 4
4700        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4701        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4702        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4703        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4704        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4705        """
4706        sh0=getShape(arg0)
4707        sh1=getShape(arg1)
4708        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4709           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4710        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4711           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4712        else:
4713            raise ValueError,"first argument must have rank 2 or 4"
4714    
4715    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4716        """
4717        generalized tensor product of transposed of arg0 and arg1:
4718    
4719        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4720    
4721        where
4722    
4723            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4724            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4725            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4726    
4727        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4728        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4729    
4730        @param arg0: first argument
4731        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4732        @param arg1: second argument
4733        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4734        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4735        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4736        """
4737        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4738        arg0,arg1=matchType(arg0,arg1)
4739        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4740        if isinstance(arg0,numarray.NumArray):
4741           if isinstance(arg1,Symbol):
4742               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4743           else:
4744               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4745                   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)
4746               arg0_c=arg0.copy()
4747               arg1_c=arg1.copy()
4748               sh0,sh1=arg0.shape,arg1.shape
4749               d0,d1,d01=1,1,1
4750               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4751               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4752               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4753               arg0_c.resize((d0,d01))
4754               arg1_c.resize((d1,d01))
4755               out=numarray.zeros((d0,d1),numarray.Float64)
4756               for i0 in range(d0):
4757                        for i1 in range(d1):
4758                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4759               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4760               return out
4761        elif isinstance(arg0,escript.Data):
4762           if isinstance(arg1,Symbol):
4763               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4764           else:
4765               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4766        else:      
4767           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4768                    
4769    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4770       """
4771       Symbol representing the general tensor product of arg0 and the transpose of arg1
4772       """
4773       def __init__(self,arg0,arg1,axis_offset=0):
4774           """
4775           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4776    
4777           @param arg0: first argument
4778           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4779           @param arg1: second argument
4780           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4781           @raise ValueError: inconsistent dimensions of arguments.
4782           @note: if both arguments have a spatial dimension, they must equal.
4783           """
4784           sh_arg0=getShape(arg0)
4785           sh_arg1=getShape(arg1)
4786           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4787           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4788           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4789           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4790           if not sh01==sh10:
4791               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)
4792           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4793    
4794       def getMyCode(self,argstrs,format="escript"):
4795          """
4796          returns a program code that can be used to evaluate the symbol.
4797    
4798          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4799          @type argstrs: C{list} of length 2 of C{str}.
4800          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4801          @type format: C{str}
4802          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4803          @rtype: C{str}
4804          @raise NotImplementedError: if the requested format is not available
4805          """
4806          if format=="escript" or format=="str" or format=="text":
4807             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4808          else:
4809             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4810    
4811       def substitute(self,argvals):
4812          """
4813          assigns new values to symbols in the definition of the symbol.
4814          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4815    
4816          @param argvals: new values assigned to symbols
4817          @type argvals: C{dict} with keywords of type L{Symbol}.
4818          @return: result of the substitution process. Operations are executed as much as possible.
4819          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4820          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4821          """
4822          if argvals.has_key(self):
4823             arg=argvals[self]
4824             if self.isAppropriateValue(arg):
4825                return arg
4826             else:
4827                raise TypeError,"%s: new value is not appropriate."%str(self)
4828          else:
4829             args=self.getSubstitutedArguments(argvals)
4830             return generalTensorTransposedProduct(args[0],args[1],args[2])
4831    
4832    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4833        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4834        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4835    
4836  #=========================================================  #=========================================================
4837  #   some little helpers  #  functions dealing with spatial dependency
4838  #=========================================================  #=========================================================
4839  def grad(arg,where=None):  def grad(arg,where=None):
4840      """      """
4841      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4842    
4843        If C{g} is the returned object, then
4844    
4845      @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.
4846                    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.
4847          - 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.
4848          - 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.
4849    
4850        @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4851        @type arg: L{escript.Data} or L{Symbol}
4852      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4853                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4854        @type where: C{None} or L{escript.FunctionSpace}
4855        @return: gradient of arg.
4856        @rtype: L{escript.Data} or L{Symbol}
4857      """      """
4858      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4859         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3617  def grad(arg,where=None): Line 4863  def grad(arg,where=None):
4863         else:         else:
4864            return arg._grad(where)            return arg._grad(where)
4865      else:      else:
4866        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4867    
4868    class Grad_Symbol(DependendSymbol):
4869       """
4870       L{Symbol} representing the result of the gradient operator
4871       """
4872       def __init__(self,arg,where=None):
4873          """
4874          initialization of gradient L{Symbol} with argument arg
4875          @param arg: argument of function
4876          @type arg: L{Symbol}.
4877          @param where: FunctionSpace in which the gradient will be calculated.
4878                      If not present or C{None} an appropriate default is used.
4879          @type where: C{None} or L{escript.FunctionSpace}
4880          """
4881          d=arg.getDim()
4882          if d==None:
4883             raise ValueError,"argument must have a spatial dimension"
4884          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4885    
4886       def getMyCode(self,argstrs,format="escript"):
4887          """
4888          returns a program code that can be used to evaluate the symbol.
4889    
4890          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4891          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4892          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4893          @type format: C{str}
4894          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4895          @rtype: C{str}
4896          @raise NotImplementedError: if the requested format is not available
4897          """
4898          if format=="escript" or format=="str"  or format=="text":
4899             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4900          else:
4901             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4902    
4903       def substitute(self,argvals):
4904          """
4905          assigns new values to symbols in the definition of the symbol.
4906          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4907    
4908          @param argvals: new values assigned to symbols
4909          @type argvals: C{dict} with keywords of type L{Symbol}.
4910          @return: result of the substitution process. Operations are executed as much as possible.
4911          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4912          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4913          """
4914          if argvals.has_key(self):
4915             arg=argvals[self]
4916             if self.isAppropriateValue(arg):
4917                return arg
4918             else:
4919                raise TypeError,"%s: new value is not appropriate."%str(self)
4920          else:
4921             arg=self.getSubstitutedArguments(argvals)
4922             return grad(arg[0],where=arg[1])
4923    
4924       def diff(self,arg):
4925          """
4926          differential of this object
4927    
4928          @param arg: the derivative is calculated with respect to arg
4929          @type arg: L{escript.Symbol}
4930          @return: derivative with respect to C{arg}
4931          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4932          """
4933          if arg==self:
4934             return identity(self.getShape())
4935          else:
4936             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4937    
4938  def integrate(arg,where=None):  def integrate(arg,where=None):
4939      """      """
4940      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}
4941      its domain.      before integration.
4942    
4943      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4944        @type arg: L{escript.Data} or L{Symbol}
4945      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4946                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4947        @type where: C{None} or L{escript.FunctionSpace}
4948        @return: integral of arg.
4949        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4950      """      """
4951      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):  
4952         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4953      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4954         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 4957  def integrate(arg,where=None):
4957         else:         else:
4958            return arg._integrate()            return arg._integrate()
4959      else:      else:
4960        raise TypeError,"integrate: Unknown argument type."         arg2=escript.Data(arg,where)
4961           if arg2.getRank()==0:
4962              return arg2._integrate()[0]
4963           else:
4964              return arg2._integrate()
4965    
4966    class Integrate_Symbol(DependendSymbol):
4967       """
4968       L{Symbol} representing the result of the spatial integration operator
4969       """
4970       def __init__(self,arg,where=None):
4971          """
4972          initialization of integration L{Symbol} with argument arg
4973          @param arg: argument of the integration
4974          @type arg: L{Symbol}.
4975          @param where: FunctionSpace in which the integration will be calculated.
4976                      If not present or C{None} an appropriate default is used.
4977          @type where: C{None} or L{escript.FunctionSpace}
4978          """
4979          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4980    
4981       def getMyCode(self,argstrs,format="escript"):
4982          """
4983          returns a program code that can be used to evaluate the symbol.
4984    
4985          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4986          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4987          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4988          @type format: C{str}
4989          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4990          @rtype: C{str}
4991          @raise NotImplementedError: if the requested format is not available
4992          """
4993          if format=="escript" or format=="str"  or format=="text":
4994             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4995          else:
4996             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4997    
4998       def substitute(self,argvals):
4999          """
5000          assigns new values to symbols in the definition of the symbol.
5001          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
5002    
5003          @param argvals: new values assigned to symbols
5004          @type argvals: C{dict} with keywords of type L{Symbol}.
5005          @return: result of the substitution process. Operations are executed as much as possible.
5006          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
5007          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
5008          """
5009          if argvals.has_key(self):
5010             arg=argvals[self]
5011             if self.isAppropriateValue(arg):
5012                return arg
5013             else:
5014                raise TypeError,"%s: new value is not appropriate."%str(self)
5015          else:
5016             arg=self.getSubstitutedArguments(argvals)
5017             return integrate(arg[0],where=arg[1])
5018    
5019       def diff(self,arg):
5020          """
5021          differential of this object
5022    
5023          @param arg: the derivative is calculated with respect to arg
5024          @type arg: L{escript.Symbol}
5025          @return: derivative with respect to C{arg}
5026          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
5027          """
5028          if arg==self:
5029             return identity(self.getShape())
5030          else:
5031             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5032    
5033    
5034  def interpolate(arg,where):  def interpolate(arg,where):
5035      """      """
5036      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where. If the argument C{arg} has the requested function space
5037        C{where} no interpolation is performed and C{arg} is returned.
5038    
5039      @param arg:    interpolant      @param arg: interpolant
5040      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
5041        @param where: FunctionSpace to be interpolated to
5042        @type where: L{escript.FunctionSpace}
5043        @return: interpolated argument
5044        @rtype: C{escript.Data} or L{Symbol}
5045      """      """
5046      if testForZero(arg):      if isinstance(arg,Symbol):
5047        return 0         return Interpolate_Symbol(arg,where)
     elif isinstance(arg,Symbol):  
        return Interpolated_Symbol(arg,where)  
5048      else:      else:
5049         return escript.Data(arg,where)         if where == arg.getFunctionSpace():
5050              return arg
5051           else:
5052              return escript.Data(arg,where)
5053    
5054  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
5055      """     """
5056      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
5057       """
5058       def __init__(self,arg,where):
5059          """
5060          initialization of interpolation L{Symbol} with argument arg
5061          @param arg: argument of the interpolation
5062          @type arg: L{Symbol}.
5063          @param where: FunctionSpace into which the argument is interpolated.
5064          @type where: L{escript.FunctionSpace}
5065          """
5066          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
5067    
5068      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
5069                    be calculated.        """
5070      @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)  
5071    
5072  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
5073      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
5074      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.
5075          @type format: C{str}
5076          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
5077          @rtype: C{str}
5078          @raise NotImplementedError: if the requested format is not available
5079          """
5080          if format=="escript" or format=="str"  or format=="text":
5081             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
5082          else:
5083             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
5084    
5085      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
5086                    to be calculated.        """
5087      """        assigns new values to symbols in the definition of the symbol.
5088      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())  
5089    
5090  #=============================        @param argvals: new values assigned to symbols
5091  #        @type argvals: C{dict} with keywords of type L{Symbol}.
5092  # 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.
5093  # 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
5094  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
5095          """
5096          if argvals.has_key(self):
5097             arg=argvals[self]
5098             if self.isAppropriateValue(arg):
5099                return arg
5100             else:
5101                raise TypeError,"%s: new value is not appropriate."%str(self)
5102          else:
5103             arg=self.getSubstitutedArguments(argvals)
5104             return interpolate(arg[0],where=arg[1])
5105    
5106  # functions involving the underlying Domain:     def diff(self,arg):
5107          """
5108          differential of this object
5109    
5110  def transpose(arg,axis=None):        @param arg: the derivative is calculated with respect to arg
5111          @type arg: L{escript.Symbol}
5112          @return: derivative with respect to C{arg}
5113          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
5114          """
5115          if arg==self:
5116             return identity(self.getShape())
5117          else:
5118             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5119    
5120    
5121    def div(arg,where=None):
5122      """      """
5123      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
5124    
5125      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
5126        @type arg: L{escript.Data} or L{Symbol}
5127        @param where: FunctionSpace in which the divergence will be calculated.
5128                      If not present or C{None} an appropriate default is used.
5129        @type where: C{None} or L{escript.FunctionSpace}
5130        @return: divergence of arg.
5131        @rtype: L{escript.Data} or L{Symbol}
5132      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
5133      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5134         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
5135      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
5136         # 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)  
5137      else:      else:
5138         return numarray.transpose(arg,axis=axis)          raise TypeError,"div: argument type not supported"
5139        if not arg.getShape()==(dim,):
5140          raise ValueError,"div: expected shape is (%s,)"%dim
5141        return trace(grad(arg,where))
5142    
5143  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
5144      """      """
5145      Return      returns the jump of arg across the continuity of the domain
5146    
5147      @param arg:      @param arg: argument
5148        @type arg: L{escript.Data} or L{Symbol}
5149        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5150                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5151        @type domain: C{None} or L{escript.Domain}
5152        @return: jump of arg
5153        @rtype: L{escript.Data} or L{Symbol}
5154      """      """
5155      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
5156         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)  
5157    
5158    def L2(arg):
5159        """
5160        returns the L2 norm of arg at where
5161        
5162        @param arg: function which L2 to be calculated.
5163        @type arg: L{escript.Data} or L{Symbol}
5164        @return: L2 norm of arg.
5165        @rtype: L{float} or L{Symbol}
5166        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5167        """
5168        return sqrt(integrate(inner(arg,arg)))
5169    #=============================
5170    #
5171    
5172  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5173      """      """
5174      resorts the component of arg according to index      resorts the component of arg according to index
5175    
5176      """      """