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

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

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

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