/[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 2100 by gross, Wed Nov 26 08:13:00 2008 UTC
# Line 1  Line 1 
1  # $Id$  
2  #  ########################################################
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
3  #  #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22  """  """
23  Utility functions for escript  Utility functions for escript
24    
 @remark:  This module is under construction and is still tested!!!  
   
25  @var __author__: name of author  @var __author__: name of author
26  @var __licence__: licence agreement  @var __copyright__: copyrights
27    @var __license__: licence agreement
28  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
29  @var __version__: version  @var __version__: version
30  @var __date__: date of the version  @var __date__: date of the version
31    @var EPSILON: smallest positive value with 1.<1.+EPSILON
32    @var DBLE_MAX: largest positive float
33  """  """
34                                                                                                                                                                                                        
35  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __licence__="contact: esys@access.uq.edu.au"  
 __url__="http://www.iservo.edu.au/esys/escript"  
 __version__="$Revision: 329 $"  
 __date__="$Date$"  
36    
37    
38  import math  import math
39  import numarray  import numarray
40  import escript  import escript
41  import os  import os
42    from esys.escript import C_GeneralTensorProduct
43  # missing tests:  from esys.escript import getVersion
44    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  
45    
46  #=========================================================  #=========================================================
47  #   some helpers:  #   some helpers:
48  #=========================================================  #=========================================================
49    def getEpsilon():
50         return escript.getMachinePrecision()
51    EPSILON=getEpsilon()
52    
53    def getMaxFloat():
54         return escript.getMaxFloat()
55    DBLE_MAX=getMaxFloat()
56    
57    
58    def getTagNames(domain):
59        """
60        returns a list of the tag names used by the domain
61    
62        
63        @param domain: a domain object
64        @type domain: L{escript.Domain}
65        @return: a list of the tag name used by the domain.
66        @rtype: C{list} of C{str}
67        """
68        return [n.strip() for n in domain.showTagNames().split(",") ]
69    
70    def insertTagNames(domain,**kwargs):
71        """
72        inserts tag names into the domain
73    
74        @param domain: a domain object
75        @type domain: C{escript.Domain}
76        @keyword <tag_name>: tag key assigned to <tag_name>
77        @type <tag_name>: C{int}
78        """
79        for  k in kwargs:
80             domain.setTagMap(k,kwargs[k])
81    
82    def insertTaggedValues(target,**kwargs):
83        """
84        inserts tagged values into the tagged using tag names
85    
86        @param target: data to be filled by tagged values
87        @type target: L{escript.Data}
88        @keyword <tag_name>: value to be used for <tag_name>
89        @type <tag_name>: C{float} or {numarray.NumArray}
90        @return: C{target}
91        @rtype: L{escript.Data}
92        """
93        for k in kwargs:
94            target.setTaggedValue(k,kwargs[k])
95        return target
96    
97  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
98      """      """
99      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.
100    
101      Example:      Example::
102    
103         tmp=Scalar(..)         tmp=Scalar(..)
104         v=Vector(..)         v=Vector(..)
105         saveVTK("solution.xml",temperature=tmp,velovity=v)         saveVTK("solution.xml",temperature=tmp,velocity=v)
106    
107      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"
108    
109      @param filename: file name of the output file      @param filename: file name of the output file
110      @type filename: C{str}      @type filename: C{str}
# Line 84  def saveVTK(filename,domain=None,**data) Line 114  def saveVTK(filename,domain=None,**data)
114      @type <name>: L{Data} object.      @type <name>: L{Data} object.
115      @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.
116      """      """
117      if domain==None:      new_data={}
118         for i in data.keys():      for n,d in data.items():
119            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
120                fs=d.getFunctionSpace()
121                domain2=fs.getDomain()
122                if fs == escript.Solution(domain2):
123                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
124                elif fs == escript.ReducedSolution(domain2):
125                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
126                else:
127                   new_data[n]=d
128                if domain==None: domain=domain2
129      if domain==None:      if domain==None:
130          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
131      else:      domain.saveVTK(filename,new_data)
         domain.saveVTK(filename,data)  
132    
133  def saveDX(filename,domain=None,**data):  def saveDX(filename,domain=None,**data):
134      """      """
135      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.
136    
137      Example:      Example::
138    
139         tmp=Scalar(..)         tmp=Scalar(..)
140         v=Vector(..)         v=Vector(..)
141         saveDX("solution.dx",temperature=tmp,velovity=v)         saveDX("solution.dx",temperature=tmp,velocity=v)
142    
143      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".
144    
145      @param filename: file name of the output file      @param filename: file name of the output file
146      @type filename: C{str}      @type filename: C{str}
# Line 112  def saveDX(filename,domain=None,**data): Line 150  def saveDX(filename,domain=None,**data):
150      @type <name>: L{Data} object.      @type <name>: L{Data} object.
151      @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.
152      """      """
153      if domain==None:      new_data={}
154         for i in data.keys():      for n,d in data.items():
155            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
156                fs=d.getFunctionSpace()
157                domain2=fs.getDomain()
158                if fs == escript.Solution(domain2):
159                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
160                elif fs == escript.ReducedSolution(domain2):
161                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
162                elif fs == escript.ContinuousFunction(domain2):
163                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
164                else:
165                   new_data[n]=d
166                if domain==None: domain=domain2
167      if domain==None:      if domain==None:
168          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
169      else:      domain.saveDX(filename,new_data)
         domain.saveDX(filename,data)  
170    
171  def kronecker(d=3):  def kronecker(d=3):
172     """     """
173     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
174    
175     @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
176     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
177     @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
178     @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}  
179     """     """
180     return identityTensor(d)     return identityTensor(d)
181    
# Line 143  def identity(shape=()): Line 190  def identity(shape=()):
190     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
191     """     """
192     if len(shape)>0:     if len(shape)>0:
193        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
194        if len(shape)==1:        if len(shape)==1:
195            for i0 in range(shape[0]):            for i0 in range(shape[0]):
196               out[i0,i0]=1.               out[i0,i0]=1.
   
197        elif len(shape)==2:        elif len(shape)==2:
198            for i0 in range(shape[0]):            for i0 in range(shape[0]):
199               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 163  def identityTensor(d=3): Line 209  def identityTensor(d=3):
209     return the dxd identity matrix     return the dxd identity matrix
210    
211     @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
212     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
213     @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
214     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
215     """     """
216     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
217        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
218     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
219           return identity((d.getDim(),))
220       else:
221           return identity((d,))
222    
223  def identityTensor4(d=3):  def identityTensor4(d=3):
224     """     """
# Line 178  def identityTensor4(d=3): Line 227  def identityTensor4(d=3):
227     @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
228     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
229     @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
230     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
231     """     """
232     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
233        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
234     return identity((d,d))     elif isinstance(d,escript.Domain):
235           return identity((d.getDim(),d.getDim()))
236       else:
237           return identity((d,d))
238    
239  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
240     """     """
# Line 191  def unitVector(i=0,d=3): Line 243  def unitVector(i=0,d=3):
243     @param i: index     @param i: index
244     @type i: C{int}     @type i: C{int}
245     @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
246     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
247     @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
248     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
249     """     """
250     return kronecker(d)[i]     return kronecker(d)[i]
251    
# Line 249  def inf(arg): Line 301  def inf(arg):
301    
302      @param arg: argument      @param arg: argument
303      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
304      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
305      @rtype: C{float}      @rtype: C{float}
306      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
307      """      """
# Line 268  def inf(arg): Line 320  def inf(arg):
320  #=========================================================================  #=========================================================================
321  #   some little helpers  #   some little helpers
322  #=========================================================================  #=========================================================================
323  def pokeShape(arg):  def getRank(arg):
324        """
325        identifies the rank of its argument
326    
327        @param arg: a given object
328        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
329        @return: the rank of the argument
330        @rtype: C{int}
331        @raise TypeError: if type of arg cannot be processed
332        """
333    
334        if isinstance(arg,numarray.NumArray):
335            return arg.rank
336        elif isinstance(arg,escript.Data):
337            return arg.getRank()
338        elif isinstance(arg,float):
339            return 0
340        elif isinstance(arg,int):
341            return 0
342        elif isinstance(arg,Symbol):
343            return arg.getRank()
344        else:
345          raise TypeError,"getShape: cannot identify shape"
346    def getShape(arg):
347      """      """
348      identifies the shape of its argument      identifies the shape of its argument
349    
# Line 290  def pokeShape(arg): Line 365  def pokeShape(arg):
365      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
366          return arg.getShape()          return arg.getShape()
367      else:      else:
368        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
369    
370  def pokeDim(arg):  def pokeDim(arg):
371      """      """
# Line 313  def commonShape(arg0,arg1): Line 388  def commonShape(arg0,arg1):
388      """      """
389      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.
390    
391      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
392      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
393      @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.
394      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
395      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
396      """      """
397      sh0=pokeShape(arg0)      sh0=getShape(arg0)
398      sh1=pokeShape(arg1)      sh1=getShape(arg1)
399      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
400         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
401               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 413  def commonDim(*args):
413      """      """
414      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.
415    
416      @param *args: given objects      @param args: given objects
417      @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
418               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
419      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 360  def testForZero(arg): Line 435  def testForZero(arg):
435    
436      @param arg: a given object      @param arg: a given object
437      @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}
438      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
439      @rtype : C{bool}      @rtype: C{bool}
440      """      """
441      try:      if isinstance(arg,numarray.NumArray):
442         return not Lsup(arg)>0.         return not Lsup(arg)>0.
443      except TypeError:      elif isinstance(arg,escript.Data):
444           return False
445        elif isinstance(arg,float):
446           return not Lsup(arg)>0.
447        elif isinstance(arg,int):
448           return not Lsup(arg)>0.
449        elif isinstance(arg,Symbol):
450           return False
451        else:
452         return False         return False
453    
454  def matchType(arg0=0.,arg1=0.):  def matchType(arg0=0.,arg1=0.):
# Line 386  def matchType(arg0=0.,arg1=0.): Line 469  def matchType(arg0=0.,arg1=0.):
469         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
470            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
471         elif isinstance(arg1,float):         elif isinstance(arg1,float):
472            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
473         elif isinstance(arg1,int):         elif isinstance(arg1,int):
474            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
475         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
476            pass            pass
477         else:         else:
# Line 412  def matchType(arg0=0.,arg1=0.): Line 495  def matchType(arg0=0.,arg1=0.):
495         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
496            pass            pass
497         elif isinstance(arg1,float):         elif isinstance(arg1,float):
498            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
499         elif isinstance(arg1,int):         elif isinstance(arg1,int):
500            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
501         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
502            pass            pass
503         else:         else:
504            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
505      elif isinstance(arg0,float):      elif isinstance(arg0,float):
506         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
507            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
508         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
509            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
510         elif isinstance(arg1,float):         elif isinstance(arg1,float):
511            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
512            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
513         elif isinstance(arg1,int):         elif isinstance(arg1,int):
514            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
515            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
516         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
517            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
518         else:         else:
519            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
520      elif isinstance(arg0,int):      elif isinstance(arg0,int):
521         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
522            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
523         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
524            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
525         elif isinstance(arg1,float):         elif isinstance(arg1,float):
526            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
527            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
528         elif isinstance(arg1,int):         elif isinstance(arg1,int):
529            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
530            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
531         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
532            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
533         else:         else:
534            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
535      else:      else:
# Line 456  def matchType(arg0=0.,arg1=0.): Line 539  def matchType(arg0=0.,arg1=0.):
539    
540  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
541      """      """
542            return representations of arg0 amd arg1 which ahve the same shape
543    
544      If shape is not given the shape "largest" shape of args is used.      @param arg0: a given object
545        @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
546      @param args: a given ob      @param arg1: a given object
547      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
548      @return: True if the argument is identical to zero.      @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
549      @rtype: C{list} of C{int}      @rtype: C{tuple}
550      """      """
551      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
552      sh0=pokeShape(arg0)      sh0=getShape(arg0)
553      sh1=pokeShape(arg1)      sh1=getShape(arg1)
554      if len(sh0)<len(sh):      if len(sh0)<len(sh):
555         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
556      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
557         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
558      else:      else:
559         return arg0,arg1         return arg0,arg1
560  #=========================================================  #=========================================================
# Line 491  class Symbol(object): Line 574  class Symbol(object):
574         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
575         symbols or any other object.         symbols or any other object.
576    
577         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
578         @type arg: C{list}         @type args: C{list}
579         @param shape: the shape         @param shape: the shape
580         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
581         @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 618  class Symbol(object):
618         """         """
619         the shape of the symbol.         the shape of the symbol.
620    
621         @return : the shape of the symbol.         @return: the shape of the symbol.
622         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
623         """         """
624         return self.__shape         return self.__shape
# Line 544  class Symbol(object): Line 627  class Symbol(object):
627         """         """
628         the spatial dimension         the spatial dimension
629    
630         @return : the spatial dimension         @return: the spatial dimension
631         @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.
632         """         """
633         return self.__dim         return self.__dim
# Line 568  class Symbol(object): Line 651  class Symbol(object):
651         """         """
652         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.
653    
654         @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].
655         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
656         @rtype: C{list} of objects         @rtype: C{list} of objects
657         @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.
658         """         """
659         out=[]         out=[]
660         for a in self.getArgument():         for a in self.getArgument():
# Line 595  class Symbol(object): Line 678  class Symbol(object):
678            if isinstance(a,Symbol):            if isinstance(a,Symbol):
679               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
680            else:            else:
681                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
682                if len(s)>0:                if len(s)>0:
683                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
684                else:                else:
685                   out.append(a)                   out.append(a)
686         return out         return out
# Line 687  class Symbol(object): Line 770  class Symbol(object):
770         else:         else:
771            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
772            if len(s)>0:            if len(s)>0:
773               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
774            else:            else:
775               return 0.               return 0.
776    
# Line 695  class Symbol(object): Line 778  class Symbol(object):
778         """         """
779         returns -self.         returns -self.
780    
781         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
782         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
783         """         """
784         return self*(-1.)         return self*(-1.)
# Line 704  class Symbol(object): Line 787  class Symbol(object):
787         """         """
788         returns +self.         returns +self.
789    
790         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
791         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
792         """         """
793         return self*(1.)         return self*(1.)
794    
795     def __abs__(self):     def __abs__(self):
796         """         """
797         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
798         """         """
799         return Abs_Symbol(self)         return Abs_Symbol(self)
800    
# Line 721  class Symbol(object): Line 804  class Symbol(object):
804    
805         @param other: object to be added to this object         @param other: object to be added to this object
806         @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}.
807         @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}
808         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
809         """         """
810         return add(self,other)         return add(self,other)
# Line 732  class Symbol(object): Line 815  class Symbol(object):
815    
816         @param other: object this object is added to         @param other: object this object is added to
817         @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}.
818         @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
819         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
820         """         """
821         return add(other,self)         return add(other,self)
# Line 743  class Symbol(object): Line 826  class Symbol(object):
826    
827         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
828         @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}.
829         @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
830         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
831         """         """
832         return add(self,-other)         return add(self,-other)
# Line 754  class Symbol(object): Line 837  class Symbol(object):
837    
838         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
839         @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}.
840         @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}.
841         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
842         """         """
843         return add(-self,other)         return add(-self,other)
# Line 765  class Symbol(object): Line 848  class Symbol(object):
848    
849         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
850         @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}.
851         @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}.
852         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
853         """         """
854         return mult(self,other)         return mult(self,other)
# Line 776  class Symbol(object): Line 859  class Symbol(object):
859    
860         @param other: object this object is multiplied with         @param other: object this object is multiplied with
861         @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}.
862         @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.
863         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
864         """         """
865         return mult(other,self)         return mult(other,self)
# Line 787  class Symbol(object): Line 870  class Symbol(object):
870    
871         @param other: object dividing this object         @param other: object dividing this object
872         @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}.
873         @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}
874         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
875         """         """
876         return quotient(self,other)         return quotient(self,other)
# Line 798  class Symbol(object): Line 881  class Symbol(object):
881    
882         @param other: object dividing this object         @param other: object dividing this object
883         @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}.
884         @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
885         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
886         """         """
887         return quotient(other,self)         return quotient(other,self)
# Line 809  class Symbol(object): Line 892  class Symbol(object):
892    
893         @param other: exponent         @param other: exponent
894         @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}.
895         @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}
896         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
897         """         """
898         return power(self,other)         return power(self,other)
# Line 820  class Symbol(object): Line 903  class Symbol(object):
903    
904         @param other: basis         @param other: basis
905         @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}.
906         @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
907         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
908         """         """
909         return power(other,self)         return power(other,self)
910    
911       def __getitem__(self,index):
912           """
913           returns the slice defined by index
914    
915           @param index: defines a
916           @type index: C{slice} or C{int} or a C{tuple} of them
917           @return: a L{Symbol} representing the slice defined by index
918           @rtype: L{DependendSymbol}
919           """
920           return GetSlice_Symbol(self,index)
921    
922  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
923     """     """
924     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.
925     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  
926        
927     Example:     Example::
928        
929     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
930     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
931     print u1==u2       print u1==u2
932     False       False
933        
934        but       but::
935    
936     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
937     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
938     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
939     print u1==u2, u1==u3       print u1==u2, u1==u3
940     True False       True False
941    
942     @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.
943     """     """
# Line 875  class DependendSymbol(Symbol): Line 969  class DependendSymbol(Symbol):
969  #=========================================================  #=========================================================
970  #  Unary operations prserving the shape  #  Unary operations prserving the shape
971  #========================================================  #========================================================
972    class GetSlice_Symbol(DependendSymbol):
973       """
974       L{Symbol} representing getting a slice for a L{Symbol}
975       """
976       def __init__(self,arg,index):
977          """
978          initialization of wherePositive L{Symbol} with argument arg
979          @param arg: argument
980          @type arg: L{Symbol}.
981          @param index: defines index
982          @type index: C{slice} or C{int} or a C{tuple} of them
983          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
984          @raises ValueError: if a step is given
985          """
986          if not isinstance(index,tuple): index=(index,)
987          if len(index)>arg.getRank():
988               raise IndexError,"GetSlice_Symbol: index out of range."
989          sh=()
990          index2=()
991          for i in range(len(index)):
992             ix=index[i]
993             if isinstance(ix,int):
994                if ix<0 or ix>=arg.getShape()[i]:
995                   raise ValueError,"GetSlice_Symbol: index out of range."
996                index2=index2+(ix,)
997             else:
998               if not ix.step==None:
999                 raise ValueError,"GetSlice_Symbol: steping is not supported."
1000               if ix.start==None:
1001                  s=0
1002               else:
1003                  s=ix.start
1004               if ix.stop==None:
1005                  e=arg.getShape()[i]
1006               else:
1007                  e=ix.stop
1008                  if e>arg.getShape()[i]:
1009                     raise IndexError,"GetSlice_Symbol: index out of range."
1010               index2=index2+(slice(s,e),)
1011               if e>s:
1012                   sh=sh+(e-s,)
1013               elif s>e:
1014                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
1015          for i in range(len(index),arg.getRank()):
1016              index2=index2+(slice(0,arg.getShape()[i]),)
1017              sh=sh+(arg.getShape()[i],)
1018          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
1019    
1020       def getMyCode(self,argstrs,format="escript"):
1021          """
1022          returns a program code that can be used to evaluate the symbol.
1023    
1024          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1025          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1026          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1027          @type format: C{str}
1028          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1029          @rtype: C{str}
1030          @raise NotImplementedError: if the requested format is not available
1031          """
1032          if format=="escript" or format=="str"  or format=="text":
1033             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
1034          else:
1035             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
1036    
1037       def substitute(self,argvals):
1038          """
1039          assigns new values to symbols in the definition of the symbol.
1040          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1041    
1042          @param argvals: new values assigned to symbols
1043          @type argvals: C{dict} with keywords of type L{Symbol}.
1044          @return: result of the substitution process. Operations are executed as much as possible.
1045          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1046          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1047          """
1048          if argvals.has_key(self):
1049             arg=argvals[self]
1050             if self.isAppropriateValue(arg):
1051                return arg
1052             else:
1053                raise TypeError,"%s: new value is not appropriate."%str(self)
1054          else:
1055             args=self.getSubstitutedArguments(argvals)
1056             arg=args[0]
1057             index=args[1]
1058             return arg.__getitem__(index)
1059    
1060  def log10(arg):  def log10(arg):
1061     """     """
1062     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
1063    
1064     @param arg: argument     @param arg: argument
1065     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1066     @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.
1067     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1068     """     """
1069     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 903  def wherePositive(arg): Line 1085  def wherePositive(arg):
1085    
1086     @param arg: argument     @param arg: argument
1087     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1088     @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.
1089     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1090     """     """
1091     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1092        if arg.rank==0:        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1093           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1094             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))  
1095     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1096        return arg._wherePositive()        return arg._wherePositive()
1097     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 953  class WherePositive_Symbol(DependendSymb Line 1131  class WherePositive_Symbol(DependendSymb
1131        @type format: C{str}        @type format: C{str}
1132        @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.
1133        @rtype: C{str}        @rtype: C{str}
1134        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1135        """        """
1136        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1137            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 989  def whereNegative(arg): Line 1167  def whereNegative(arg):
1167    
1168     @param arg: argument     @param arg: argument
1169     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1170     @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.
1171     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1172     """     """
1173     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1174        if arg.rank==0:        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1175           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1176             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))  
1177     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1178        return arg._whereNegative()        return arg._whereNegative()
1179     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1039  class WhereNegative_Symbol(DependendSymb Line 1213  class WhereNegative_Symbol(DependendSymb
1213        @type format: C{str}        @type format: C{str}
1214        @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.
1215        @rtype: C{str}        @rtype: C{str}
1216        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1217        """        """
1218        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1219            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1075  def whereNonNegative(arg): Line 1249  def whereNonNegative(arg):
1249    
1250     @param arg: argument     @param arg: argument
1251     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1252     @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.
1253     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1254     """     """
1255     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1256        if arg.rank==0:        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1257           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1258             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))  
1259     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1260        return arg._whereNonNegative()        return arg._whereNonNegative()
1261     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1109  def whereNonPositive(arg): Line 1279  def whereNonPositive(arg):
1279    
1280     @param arg: argument     @param arg: argument
1281     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1282     @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.
1283     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1284     """     """
1285     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1286        if arg.rank==0:        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1287           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1288             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.  
1289     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1290        return arg._whereNonPositive()        return arg._whereNonPositive()
1291     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1145  def whereZero(arg,tol=0.): Line 1311  def whereZero(arg,tol=0.):
1311     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1312     @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.
1313     @type tol: C{float}     @type tol: C{float}
1314     @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.
1315     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1316     """     """
1317     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1318        if arg.rank==0:        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1319           if abs(arg)<=tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1320             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.  
1321     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1322        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1323     elif isinstance(arg,float):     elif isinstance(arg,float):
1324        if abs(arg)<=tol:        if abs(arg)<=tol:
1325          return 1.          return 1.
# Line 1198  class WhereZero_Symbol(DependendSymbol): Line 1357  class WhereZero_Symbol(DependendSymbol):
1357        @type format: C{str}        @type format: C{str}
1358        @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.
1359        @rtype: C{str}        @rtype: C{str}
1360        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1361        """        """
1362        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1363           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 1391  def whereNonZero(arg,tol=0.):
1391    
1392     @param arg: argument     @param arg: argument
1393     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1394     @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.
1395     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1396     """     """
1397     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1398        if arg.rank==0:        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1399          if abs(arg)>tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1400             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.  
1401     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1402        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1403     elif isinstance(arg,float):     elif isinstance(arg,float):
1404        if abs(arg)>tol:        if abs(arg)>tol:
1405          return 1.          return 1.
# Line 1263  def whereNonZero(arg,tol=0.): Line 1415  def whereNonZero(arg,tol=0.):
1415     else:     else:
1416        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1417    
1418    def erf(arg):
1419       """
1420       returns erf of argument arg
1421    
1422       @param arg: argument
1423       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1424       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1425       @raises TypeError: if the type of the argument is not expected.
1426       """
1427       if isinstance(arg,escript.Data):
1428          return arg._erf()
1429       else:
1430          raise TypeError,"erf: Unknown argument type."
1431    
1432  def sin(arg):  def sin(arg):
1433     """     """
1434     returns sine of argument arg     returns sine of argument arg
1435    
1436     @param arg: argument     @param arg: argument
1437     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1438     @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.
1439     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1440     """     """
1441     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1307  class Sin_Symbol(DependendSymbol): Line 1473  class Sin_Symbol(DependendSymbol):
1473        @type format: C{str}        @type format: C{str}
1474        @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.
1475        @rtype: C{str}        @rtype: C{str}
1476        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1477        """        """
1478        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1479            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1359  def cos(arg): Line 1525  def cos(arg):
1525    
1526     @param arg: argument     @param arg: argument
1527     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1528     @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.
1529     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1530     """     """
1531     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1397  class Cos_Symbol(DependendSymbol): Line 1563  class Cos_Symbol(DependendSymbol):
1563        @type format: C{str}        @type format: C{str}
1564        @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.
1565        @rtype: C{str}        @rtype: C{str}
1566        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1567        """        """
1568        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1569            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1449  def tan(arg): Line 1615  def tan(arg):
1615    
1616     @param arg: argument     @param arg: argument
1617     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1618     @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.
1619     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1620     """     """
1621     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1487  class Tan_Symbol(DependendSymbol): Line 1653  class Tan_Symbol(DependendSymbol):
1653        @type format: C{str}        @type format: C{str}
1654        @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.
1655        @rtype: C{str}        @rtype: C{str}
1656        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1657        """        """
1658        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1659            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1539  def asin(arg): Line 1705  def asin(arg):
1705    
1706     @param arg: argument     @param arg: argument
1707     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1708     @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.
1709     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1710     """     """
1711     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1577  class Asin_Symbol(DependendSymbol): Line 1743  class Asin_Symbol(DependendSymbol):
1743        @type format: C{str}        @type format: C{str}
1744        @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.
1745        @rtype: C{str}        @rtype: C{str}
1746        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1747        """        """
1748        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1749            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1629  def acos(arg): Line 1795  def acos(arg):
1795    
1796     @param arg: argument     @param arg: argument
1797     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1798     @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.
1799     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1800     """     """
1801     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1667  class Acos_Symbol(DependendSymbol): Line 1833  class Acos_Symbol(DependendSymbol):
1833        @type format: C{str}        @type format: C{str}
1834        @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.
1835        @rtype: C{str}        @rtype: C{str}
1836        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1837        """        """
1838        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1839            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1719  def atan(arg): Line 1885  def atan(arg):
1885    
1886     @param arg: argument     @param arg: argument
1887     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1888     @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.
1889     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1890     """     """
1891     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1757  class Atan_Symbol(DependendSymbol): Line 1923  class Atan_Symbol(DependendSymbol):
1923        @type format: C{str}        @type format: C{str}
1924        @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.
1925        @rtype: C{str}        @rtype: C{str}
1926        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1927        """        """
1928        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1929            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1809  def sinh(arg): Line 1975  def sinh(arg):
1975    
1976     @param arg: argument     @param arg: argument
1977     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1978     @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.
1979     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1980     """     """
1981     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1847  class Sinh_Symbol(DependendSymbol): Line 2013  class Sinh_Symbol(DependendSymbol):
2013        @type format: C{str}        @type format: C{str}
2014        @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.
2015        @rtype: C{str}        @rtype: C{str}
2016        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2017        """        """
2018        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2019            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1899  def cosh(arg): Line 2065  def cosh(arg):
2065    
2066     @param arg: argument     @param arg: argument
2067     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2068     @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.
2069     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2070     """     """
2071     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1937  class Cosh_Symbol(DependendSymbol): Line 2103  class Cosh_Symbol(DependendSymbol):
2103        @type format: C{str}        @type format: C{str}
2104        @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.
2105        @rtype: C{str}        @rtype: C{str}
2106        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2107        """        """
2108        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2109            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1989  def tanh(arg): Line 2155  def tanh(arg):
2155    
2156     @param arg: argument     @param arg: argument
2157     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2158     @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.
2159     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2160     """     """
2161     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2027  class Tanh_Symbol(DependendSymbol): Line 2193  class Tanh_Symbol(DependendSymbol):
2193        @type format: C{str}        @type format: C{str}
2194        @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.
2195        @rtype: C{str}        @rtype: C{str}
2196        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2197        """        """
2198        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2199            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2079  def asinh(arg): Line 2245  def asinh(arg):
2245    
2246     @param arg: argument     @param arg: argument
2247     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2248     @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.
2249     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2250     """     """
2251     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2117  class Asinh_Symbol(DependendSymbol): Line 2283  class Asinh_Symbol(DependendSymbol):
2283        @type format: C{str}        @type format: C{str}
2284        @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.
2285        @rtype: C{str}        @rtype: C{str}
2286        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2287        """        """
2288        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2289            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2169  def acosh(arg): Line 2335  def acosh(arg):
2335    
2336     @param arg: argument     @param arg: argument
2337     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2338     @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.
2339     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2340     """     """
2341     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2207  class Acosh_Symbol(DependendSymbol): Line 2373  class Acosh_Symbol(DependendSymbol):
2373        @type format: C{str}        @type format: C{str}
2374        @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.
2375        @rtype: C{str}        @rtype: C{str}
2376        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2377        """        """
2378        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2379            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2259  def atanh(arg): Line 2425  def atanh(arg):
2425    
2426     @param arg: argument     @param arg: argument
2427     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2428     @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.
2429     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2430     """     """
2431     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2297  class Atanh_Symbol(DependendSymbol): Line 2463  class Atanh_Symbol(DependendSymbol):
2463        @type format: C{str}        @type format: C{str}
2464        @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.
2465        @rtype: C{str}        @rtype: C{str}
2466        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2467        """        """
2468        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2469            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2349  def exp(arg): Line 2515  def exp(arg):
2515    
2516     @param arg: argument     @param arg: argument
2517     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2518     @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.
2519     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2520     """     """
2521     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2387  class Exp_Symbol(DependendSymbol): Line 2553  class Exp_Symbol(DependendSymbol):
2553        @type format: C{str}        @type format: C{str}
2554        @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.
2555        @rtype: C{str}        @rtype: C{str}
2556        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2557        """        """
2558        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2559            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2439  def sqrt(arg): Line 2605  def sqrt(arg):
2605    
2606     @param arg: argument     @param arg: argument
2607     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2608     @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.
2609     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2610     """     """
2611     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2477  class Sqrt_Symbol(DependendSymbol): Line 2643  class Sqrt_Symbol(DependendSymbol):
2643        @type format: C{str}        @type format: C{str}
2644        @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.
2645        @rtype: C{str}        @rtype: C{str}
2646        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2647        """        """
2648        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2649            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2529  def log(arg): Line 2695  def log(arg):
2695    
2696     @param arg: argument     @param arg: argument
2697     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2698     @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.
2699     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2700     """     """
2701     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2567  class Log_Symbol(DependendSymbol): Line 2733  class Log_Symbol(DependendSymbol):
2733        @type format: C{str}        @type format: C{str}
2734        @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.
2735        @rtype: C{str}        @rtype: C{str}
2736        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2737        """        """
2738        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2739            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2619  def sign(arg): Line 2785  def sign(arg):
2785    
2786     @param arg: argument     @param arg: argument
2787     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2788     @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.
2789     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2790     """     """
2791     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2667  class Abs_Symbol(DependendSymbol): Line 2833  class Abs_Symbol(DependendSymbol):
2833        @type format: C{str}        @type format: C{str}
2834        @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.
2835        @rtype: C{str}        @rtype: C{str}
2836        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2837        """        """
2838        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2839            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2719  def minval(arg): Line 2885  def minval(arg):
2885    
2886     @param arg: argument     @param arg: argument
2887     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2888     @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.
2889     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2890     """     """
2891     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2760  class Minval_Symbol(DependendSymbol): Line 2926  class Minval_Symbol(DependendSymbol):
2926        @type format: C{str}        @type format: C{str}
2927        @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.
2928        @rtype: C{str}        @rtype: C{str}
2929        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2930        """        """
2931        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2932            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2796  def maxval(arg): Line 2962  def maxval(arg):
2962    
2963     @param arg: argument     @param arg: argument
2964     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2965     @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.
2966     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2967     """     """
2968     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2837  class Maxval_Symbol(DependendSymbol): Line 3003  class Maxval_Symbol(DependendSymbol):
3003        @type format: C{str}        @type format: C{str}
3004        @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.
3005        @rtype: C{str}        @rtype: C{str}
3006        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3007        """        """
3008        if isinstance(argstrs,list):        if isinstance(argstrs,list):
3009            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2873  def length(arg): Line 3039  def length(arg):
3039    
3040     @param arg: argument     @param arg: argument
3041     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3042     @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.
3043     """     """
3044     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
3045    
3046    def trace(arg,axis_offset=0):
3047       """
3048       returns the trace of arg which the sum of arg[k,k] over k.
3049    
3050       @param arg: argument
3051       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3052       @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
3053                      C{axis_offset} and axis_offset+1 must be equal.
3054       @type axis_offset: C{int}
3055       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
3056       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3057       """
3058       if isinstance(arg,numarray.NumArray):
3059          sh=arg.shape
3060          if len(sh)<2:
3061            raise ValueError,"rank of argument must be greater than 1"
3062          if axis_offset<0 or axis_offset>len(sh)-2:
3063            raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
3064          s1=1
3065          for i in range(axis_offset): s1*=sh[i]
3066          s2=1
3067          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3068          if not sh[axis_offset] == sh[axis_offset+1]:
3069            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3070          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
3071          out=numarray.zeros([s1,s2],numarray.Float64)
3072          for i1 in range(s1):
3073            for i2 in range(s2):
3074                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
3075          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3076          return out
3077       elif isinstance(arg,escript.Data):
3078          if arg.getRank()<2:
3079            raise ValueError,"rank of argument must be greater than 1"
3080          if axis_offset<0 or axis_offset>arg.getRank()-2:
3081            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3082          s=list(arg.getShape())        
3083          if not s[axis_offset] == s[axis_offset+1]:
3084            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3085          return arg._trace(axis_offset)
3086       elif isinstance(arg,float):
3087          raise TypeError,"illegal argument type float."
3088       elif isinstance(arg,int):
3089          raise TypeError,"illegal argument type int."
3090       elif isinstance(arg,Symbol):
3091          return Trace_Symbol(arg,axis_offset)
3092       else:
3093          raise TypeError,"Unknown argument type."
3094    
3095    class Trace_Symbol(DependendSymbol):
3096       """
3097       L{Symbol} representing the result of the trace function
3098       """
3099       def __init__(self,arg,axis_offset=0):
3100          """
3101          initialization of trace L{Symbol} with argument arg
3102          @param arg: argument of function
3103          @type arg: L{Symbol}.
3104          @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
3105                      C{axis_offset} and axis_offset+1 must be equal.
3106          @type axis_offset: C{int}
3107          """
3108          if arg.getRank()<2:
3109            raise ValueError,"rank of argument must be greater than 1"
3110          if axis_offset<0 or axis_offset>arg.getRank()-2:
3111            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3112          s=list(arg.getShape())        
3113          if not s[axis_offset] == s[axis_offset+1]:
3114            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3115          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3116    
3117       def getMyCode(self,argstrs,format="escript"):
3118          """
3119          returns a program code that can be used to evaluate the symbol.
3120    
3121          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3122          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3123          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3124          @type format: C{str}
3125          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3126          @rtype: C{str}
3127          @raise NotImplementedError: if the requested format is not available
3128          """
3129          if format=="escript" or format=="str"  or format=="text":
3130             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3131          else:
3132             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3133    
3134       def substitute(self,argvals):
3135          """
3136          assigns new values to symbols in the definition of the symbol.
3137          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3138    
3139          @param argvals: new values assigned to symbols
3140          @type argvals: C{dict} with keywords of type L{Symbol}.
3141          @return: result of the substitution process. Operations are executed as much as possible.
3142          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3143          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3144          """
3145          if argvals.has_key(self):
3146             arg=argvals[self]
3147             if self.isAppropriateValue(arg):
3148                return arg
3149             else:
3150                raise TypeError,"%s: new value is not appropriate."%str(self)
3151          else:
3152             arg=self.getSubstitutedArguments(argvals)
3153             return trace(arg[0],axis_offset=arg[1])
3154    
3155       def diff(self,arg):
3156          """
3157          differential of this object
3158    
3159          @param arg: the derivative is calculated with respect to arg
3160          @type arg: L{escript.Symbol}
3161          @return: derivative with respect to C{arg}
3162          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3163          """
3164          if arg==self:
3165             return identity(self.getShape())
3166          else:
3167             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3168    
3169    def transpose(arg,axis_offset=None):
3170       """
3171       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3172    
3173       @param arg: argument
3174       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3175       @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.
3176                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3177       @type axis_offset: C{int}
3178       @return: transpose of arg
3179       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3180       """
3181       if isinstance(arg,numarray.NumArray):
3182          if axis_offset==None: axis_offset=int(arg.rank/2)
3183          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3184       elif isinstance(arg,escript.Data):
3185          r=arg.getRank()
3186          if axis_offset==None: axis_offset=int(r/2)
3187          if axis_offset<0 or axis_offset>r:
3188            raise ValueError,"axis_offset must be between 0 and %s"%r
3189          return arg._transpose(axis_offset)
3190       elif isinstance(arg,float):
3191          if not ( axis_offset==0 or axis_offset==None):
3192            raise ValueError,"axis_offset must be 0 for float argument"
3193          return arg
3194       elif isinstance(arg,int):
3195          if not ( axis_offset==0 or axis_offset==None):
3196            raise ValueError,"axis_offset must be 0 for int argument"
3197          return float(arg)
3198       elif isinstance(arg,Symbol):
3199          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3200          return Transpose_Symbol(arg,axis_offset)
3201       else:
3202          raise TypeError,"Unknown argument type."
3203    
3204    class Transpose_Symbol(DependendSymbol):
3205       """
3206       L{Symbol} representing the result of the transpose function
3207       """
3208       def __init__(self,arg,axis_offset=None):
3209          """
3210          initialization of transpose L{Symbol} with argument arg
3211    
3212          @param arg: argument of function
3213          @type arg: L{Symbol}.
3214          @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.
3215                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3216          @type axis_offset: C{int}
3217          """
3218          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3219          if axis_offset<0 or axis_offset>arg.getRank():
3220            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3221          s=arg.getShape()
3222          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3223    
3224       def getMyCode(self,argstrs,format="escript"):
3225          """
3226          returns a program code that can be used to evaluate the symbol.
3227    
3228          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3229          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3230          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3231          @type format: C{str}
3232          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3233          @rtype: C{str}
3234          @raise NotImplementedError: if the requested format is not available
3235          """
3236          if format=="escript" or format=="str"  or format=="text":
3237             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3238          else:
3239             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3240    
3241       def substitute(self,argvals):
3242          """
3243          assigns new values to symbols in the definition of the symbol.
3244          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3245    
3246          @param argvals: new values assigned to symbols
3247          @type argvals: C{dict} with keywords of type L{Symbol}.
3248          @return: result of the substitution process. Operations are executed as much as possible.
3249          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3250          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3251          """
3252          if argvals.has_key(self):
3253             arg=argvals[self]
3254             if self.isAppropriateValue(arg):
3255                return arg
3256             else:
3257                raise TypeError,"%s: new value is not appropriate."%str(self)
3258          else:
3259             arg=self.getSubstitutedArguments(argvals)
3260             return transpose(arg[0],axis_offset=arg[1])
3261    
3262       def diff(self,arg):
3263          """
3264          differential of this object
3265    
3266          @param arg: the derivative is calculated with respect to arg
3267          @type arg: L{escript.Symbol}
3268          @return: derivative with respect to C{arg}
3269          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3270          """
3271          if arg==self:
3272             return identity(self.getShape())
3273          else:
3274             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3275    
3276    def swap_axes(arg,axis0=0,axis1=1):
3277       """
3278       returns the swap of arg by swaping the components axis0 and axis1
3279    
3280       @param arg: argument
3281       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3282       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3283       @type axis0: C{int}
3284       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3285       @type axis1: C{int}
3286       @return: C{arg} with swaped components
3287       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3288       """
3289       if axis0 > axis1:
3290          axis0,axis1=axis1,axis0
3291       if isinstance(arg,numarray.NumArray):
3292          return numarray.swapaxes(arg,axis0,axis1)
3293       elif isinstance(arg,escript.Data):
3294          return arg._swap_axes(axis0,axis1)
3295       elif isinstance(arg,float):
3296          raise TyepError,"float argument is not supported."
3297       elif isinstance(arg,int):
3298          raise TyepError,"int argument is not supported."
3299       elif isinstance(arg,Symbol):
3300          return SwapAxes_Symbol(arg,axis0,axis1)
3301       else:
3302          raise TypeError,"Unknown argument type."
3303    
3304    class SwapAxes_Symbol(DependendSymbol):
3305       """
3306       L{Symbol} representing the result of the swap function
3307       """
3308       def __init__(self,arg,axis0=0,axis1=1):
3309          """
3310          initialization of swap L{Symbol} with argument arg
3311    
3312          @param arg: argument
3313          @type arg: L{Symbol}.
3314          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3315          @type axis0: C{int}
3316          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3317          @type axis1: C{int}
3318          """
3319          if arg.getRank()<2:
3320             raise ValueError,"argument must have at least rank 2."
3321          if axis0<0 or axis0>arg.getRank()-1:
3322             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3323          if axis1<0 or axis1>arg.getRank()-1:
3324             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3325          if axis0 == axis1:
3326             raise ValueError,"axis indices must be different."
3327          if axis0 > axis1:
3328             axis0,axis1=axis1,axis0
3329          s=arg.getShape()
3330          s_out=[]
3331          for i in range(len(s)):
3332             if i == axis0:
3333                s_out.append(s[axis1])
3334             elif i == axis1:
3335                s_out.append(s[axis0])
3336             else:
3337                s_out.append(s[i])
3338          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3339    
3340       def getMyCode(self,argstrs,format="escript"):
3341          """
3342          returns a program code that can be used to evaluate the symbol.
3343    
3344          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3345          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3346          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3347          @type format: C{str}
3348          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3349          @rtype: C{str}
3350          @raise NotImplementedError: if the requested format is not available
3351          """
3352          if format=="escript" or format=="str"  or format=="text":
3353             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3354          else:
3355             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3356    
3357       def substitute(self,argvals):
3358          """
3359          assigns new values to symbols in the definition of the symbol.
3360          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3361    
3362          @param argvals: new values assigned to symbols
3363          @type argvals: C{dict} with keywords of type L{Symbol}.
3364          @return: result of the substitution process. Operations are executed as much as possible.
3365          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3366          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3367          """
3368          if argvals.has_key(self):
3369             arg=argvals[self]
3370             if self.isAppropriateValue(arg):
3371                return arg
3372             else:
3373                raise TypeError,"%s: new value is not appropriate."%str(self)
3374          else:
3375             arg=self.getSubstitutedArguments(argvals)
3376             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3377    
3378       def diff(self,arg):
3379          """
3380          differential of this object
3381    
3382          @param arg: the derivative is calculated with respect to arg
3383          @type arg: L{escript.Symbol}
3384          @return: derivative with respect to C{arg}
3385          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3386          """
3387          if arg==self:
3388             return identity(self.getShape())
3389          else:
3390             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3391    
3392    def symmetric(arg):
3393        """
3394        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3395    
3396        @param arg: square matrix. Must have rank 2 or 4 and be square.
3397        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3398        @return: symmetric part of arg
3399        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3400        """
3401        if isinstance(arg,numarray.NumArray):
3402          if arg.rank==2:
3403            if not (arg.shape[0]==arg.shape[1]):
3404               raise ValueError,"argument must be square."
3405          elif arg.rank==4:
3406            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3407               raise ValueError,"argument must be square."
3408          else:
3409            raise ValueError,"rank 2 or 4 is required."
3410          return (arg+transpose(arg))/2
3411        elif isinstance(arg,escript.Data):
3412          if arg.getRank()==2:
3413            if not (arg.getShape()[0]==arg.getShape()[1]):
3414               raise ValueError,"argument must be square."
3415            return arg._symmetric()
3416          elif arg.getRank()==4:
3417            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3418               raise ValueError,"argument must be square."
3419            return arg._symmetric()
3420          else:
3421            raise ValueError,"rank 2 or 4 is required."
3422        elif isinstance(arg,float):
3423          return arg
3424        elif isinstance(arg,int):
3425          return float(arg)
3426        elif isinstance(arg,Symbol):
3427          if arg.getRank()==2:
3428            if not (arg.getShape()[0]==arg.getShape()[1]):
3429               raise ValueError,"argument must be square."
3430          elif arg.getRank()==4:
3431            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3432               raise ValueError,"argument must be square."
3433          else:
3434            raise ValueError,"rank 2 or 4 is required."
3435          return (arg+transpose(arg))/2
3436        else:
3437          raise TypeError,"symmetric: Unknown argument type."
3438    
3439    def nonsymmetric(arg):
3440        """
3441        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3442    
3443        @param arg: square matrix. Must have rank 2 or 4 and be square.
3444        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3445        @return: nonsymmetric part of arg
3446        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3447        """
3448        if isinstance(arg,numarray.NumArray):
3449          if arg.rank==2:
3450            if not (arg.shape[0]==arg.shape[1]):
3451               raise ValueError,"nonsymmetric: argument must be square."
3452          elif arg.rank==4:
3453            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3454               raise ValueError,"nonsymmetric: argument must be square."
3455          else:
3456            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3457          return (arg-transpose(arg))/2
3458        elif isinstance(arg,escript.Data):
3459          if arg.getRank()==2:
3460            if not (arg.getShape()[0]==arg.getShape()[1]):
3461               raise ValueError,"argument must be square."
3462            return arg._nonsymmetric()
3463          elif arg.getRank()==4:
3464            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3465               raise ValueError,"argument must be square."
3466            return arg._nonsymmetric()
3467          else:
3468            raise ValueError,"rank 2 or 4 is required."
3469        elif isinstance(arg,float):
3470          return arg
3471        elif isinstance(arg,int):
3472          return float(arg)
3473        elif isinstance(arg,Symbol):
3474          if arg.getRank()==2:
3475            if not (arg.getShape()[0]==arg.getShape()[1]):
3476               raise ValueError,"nonsymmetric: argument must be square."
3477          elif arg.getRank()==4:
3478            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3479               raise ValueError,"nonsymmetric: argument must be square."
3480          else:
3481            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3482          return (arg-transpose(arg))/2
3483        else:
3484          raise TypeError,"nonsymmetric: Unknown argument type."
3485    
3486    def inverse(arg):
3487        """
3488        returns the inverse of the square matrix arg.
3489    
3490        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3491        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3492        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3493        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3494        @note: for L{escript.Data} objects the dimension is restricted to 3.
3495        """
3496        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3497        if isinstance(arg,numarray.NumArray):
3498          return numarray.linear_algebra.inverse(arg)
3499        elif isinstance(arg,escript.Data):
3500          return escript_inverse(arg)
3501        elif isinstance(arg,float):
3502          return 1./arg
3503        elif isinstance(arg,int):
3504          return 1./float(arg)
3505        elif isinstance(arg,Symbol):
3506          return Inverse_Symbol(arg)
3507        else:
3508          raise TypeError,"inverse: Unknown argument type."
3509    
3510    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3511          "arg is a Data objects!!!"
3512          if not arg.getRank()==2:
3513            raise ValueError,"escript_inverse: argument must have rank 2"
3514          s=arg.getShape()      
3515          if not s[0] == s[1]:
3516            raise ValueError,"escript_inverse: argument must be a square matrix."
3517          out=escript.Data(0.,s,arg.getFunctionSpace())
3518          if s[0]==1:
3519              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3520                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3521              out[0,0]=1./arg[0,0]
3522          elif s[0]==2:
3523              A11=arg[0,0]
3524              A12=arg[0,1]
3525              A21=arg[1,0]
3526              A22=arg[1,1]
3527              D = A11*A22-A12*A21
3528              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3529                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3530              D=1./D
3531              out[0,0]= A22*D
3532              out[1,0]=-A21*D
3533              out[0,1]=-A12*D
3534              out[1,1]= A11*D
3535          elif s[0]==3:
3536              A11=arg[0,0]
3537              A21=arg[1,0]
3538              A31=arg[2,0]
3539              A12=arg[0,1]
3540              A22=arg[1,1]
3541              A32=arg[2,1]
3542              A13=arg[0,2]
3543              A23=arg[1,2]
3544              A33=arg[2,2]
3545              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3546              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3547                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3548              D=1./D
3549              out[0,0]=(A22*A33-A23*A32)*D
3550              out[1,0]=(A31*A23-A21*A33)*D
3551              out[2,0]=(A21*A32-A31*A22)*D
3552              out[0,1]=(A13*A32-A12*A33)*D
3553              out[1,1]=(A11*A33-A31*A13)*D
3554              out[2,1]=(A12*A31-A11*A32)*D
3555              out[0,2]=(A12*A23-A13*A22)*D
3556              out[1,2]=(A13*A21-A11*A23)*D
3557              out[2,2]=(A11*A22-A12*A21)*D
3558          else:
3559             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3560          return out
3561    
3562    class Inverse_Symbol(DependendSymbol):
3563       """
3564       L{Symbol} representing the result of the inverse function
3565       """
3566       def __init__(self,arg):
3567          """
3568          initialization of inverse L{Symbol} with argument arg
3569          @param arg: argument of function
3570          @type arg: L{Symbol}.
3571          """
3572          if not arg.getRank()==2:
3573            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3574          s=arg.getShape()
3575          if not s[0] == s[1]:
3576            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3577          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3578    
3579       def getMyCode(self,argstrs,format="escript"):
3580          """
3581          returns a program code that can be used to evaluate the symbol.
3582    
3583          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3584          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3585          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3586          @type format: C{str}
3587          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3588          @rtype: C{str}
3589          @raise NotImplementedError: if the requested format is not available
3590          """
3591          if format=="escript" or format=="str"  or format=="text":
3592             return "inverse(%s)"%argstrs[0]
3593          else:
3594             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3595    
3596       def substitute(self,argvals):
3597          """
3598          assigns new values to symbols in the definition of the symbol.
3599          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3600    
3601          @param argvals: new values assigned to symbols
3602          @type argvals: C{dict} with keywords of type L{Symbol}.
3603          @return: result of the substitution process. Operations are executed as much as possible.
3604          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3605          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3606          """
3607          if argvals.has_key(self):
3608             arg=argvals[self]
3609             if self.isAppropriateValue(arg):
3610                return arg
3611             else:
3612                raise TypeError,"%s: new value is not appropriate."%str(self)
3613          else:
3614             arg=self.getSubstitutedArguments(argvals)
3615             return inverse(arg[0])
3616    
3617       def diff(self,arg):
3618          """
3619          differential of this object
3620    
3621          @param arg: the derivative is calculated with respect to arg
3622          @type arg: L{escript.Symbol}
3623          @return: derivative with respect to C{arg}
3624          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3625          """
3626          if arg==self:
3627             return identity(self.getShape())
3628          else:
3629             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3630    
3631    def eigenvalues(arg):
3632        """
3633        returns the eigenvalues of the square matrix arg.
3634    
3635        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3636                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3637        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3638        @return: the eigenvalues in increasing order.
3639        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3640        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3641        """
3642        if isinstance(arg,numarray.NumArray):
3643          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3644          out.sort()
3645          return out
3646        elif isinstance(arg,escript.Data):
3647          return arg._eigenvalues()
3648        elif isinstance(arg,Symbol):
3649          if not arg.getRank()==2:
3650            raise ValueError,"eigenvalues: argument must have rank 2"
3651          s=arg.getShape()      
3652          if not s[0] == s[1]:
3653            raise ValueError,"eigenvalues: argument must be a square matrix."
3654          if s[0]==1:
3655              return arg[0]
3656          elif s[0]==2:
3657              arg1=symmetric(arg)
3658              A11=arg1[0,0]
3659              A12=arg1[0,1]
3660              A22=arg1[1,1]
3661              trA=(A11+A22)/2.
3662              A11-=trA
3663              A22-=trA
3664              s=sqrt(A12**2-A11*A22)
3665              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3666          elif s[0]==3:
3667              arg1=symmetric(arg)
3668              A11=arg1[0,0]
3669              A12=arg1[0,1]
3670              A22=arg1[1,1]
3671              A13=arg1[0,2]
3672              A23=arg1[1,2]
3673              A33=arg1[2,2]
3674              trA=(A11+A22+A33)/3.
3675              A11-=trA
3676              A22-=trA
3677              A33-=trA
3678              A13_2=A13**2
3679              A23_2=A23**2
3680              A12_2=A12**2
3681              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3682              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3683              sq_p=sqrt(p/3.)
3684              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
3685              sq_p*=2.
3686              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3687               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3688               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3689              return trA+sq_p*f
3690          else:
3691             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3692        elif isinstance(arg,float):
3693          return arg
3694        elif isinstance(arg,int):
3695          return float(arg)
3696        else:
3697          raise TypeError,"eigenvalues: Unknown argument type."
3698    
3699    def eigenvalues_and_eigenvectors(arg):
3700        """
3701        returns the eigenvalues and eigenvectors of the square matrix arg.
3702    
3703        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3704                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3705        @type arg: L{escript.Data}
3706        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3707                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3708                 the eigenvector coresponding to the i-th eigenvalue.
3709        @rtype: L{tuple} of L{escript.Data}.
3710        @note: The dimension is restricted to 3.
3711        """
3712        if isinstance(arg,numarray.NumArray):
3713          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3714        elif isinstance(arg,escript.Data):
3715          return arg._eigenvalues_and_eigenvectors()
3716        elif isinstance(arg,Symbol):
3717          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3718        elif isinstance(arg,float):
3719          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3720        elif isinstance(arg,int):
3721          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3722        else:
3723          raise TypeError,"eigenvalues: Unknown argument type."
3724  #=======================================================  #=======================================================
3725  #  Binary operations:  #  Binary operations:
3726  #=======================================================  #=======================================================
# Line 2920  class Add_Symbol(DependendSymbol): Line 3764  class Add_Symbol(DependendSymbol):
3764         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3765         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3766         """         """
3767         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3768         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3769         if not sh0==sh1:         if not sh0==sh1:
3770            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3771         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 3780  class Add_Symbol(DependendSymbol):
3780        @type format: C{str}        @type format: C{str}
3781        @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.
3782        @rtype: C{str}        @rtype: C{str}
3783        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3784        """        """
3785        if format=="str" or format=="text":        if format=="str" or format=="text":
3786           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 2964  class Add_Symbol(DependendSymbol): Line 3808  class Add_Symbol(DependendSymbol):
3808              raise TypeError,"%s: new value is not appropriate."%str(self)              raise TypeError,"%s: new value is not appropriate."%str(self)
3809        else:        else:
3810           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
3811           return add(args[0],args[1])           out=add(args[0],args[1])
3812             return out
3813    
3814     def diff(self,arg):     def diff(self,arg):
3815        """        """
# Line 2995  def mult(arg0,arg1): Line 3840  def mult(arg0,arg1):
3840         """         """
3841         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3842         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3843            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3844         else:         else:
3845            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3846                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3019  class Mult_Symbol(DependendSymbol): Line 3864  class Mult_Symbol(DependendSymbol):
3864         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3865         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3866         """         """
3867         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3868         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3869         if not sh0==sh1:         if not sh0==sh1:
3870            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3871         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 3880  class Mult_Symbol(DependendSymbol):
3880        @type format: C{str}        @type format: C{str}
3881        @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.
3882        @rtype: C{str}        @rtype: C{str}
3883        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3884        """        """
3885        if format=="str" or format=="text":        if format=="str" or format=="text":
3886           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3095  def quotient(arg0,arg1): Line 3940  def quotient(arg0,arg1):
3940         """         """
3941         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3942         if testForZero(args[0]):         if testForZero(args[0]):
3943            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3944         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3945            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3946               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3124  class Quotient_Symbol(DependendSymbol): Line 3969  class Quotient_Symbol(DependendSymbol):
3969         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3970         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3971         """         """
3972         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3973         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3974         if not sh0==sh1:         if not sh0==sh1:
3975            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3976         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 3985  class Quotient_Symbol(DependendSymbol):
3985        @type format: C{str}        @type format: C{str}
3986        @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.
3987        @rtype: C{str}        @rtype: C{str}
3988        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3989        """        """
3990        if format=="str" or format=="text":        if format=="str" or format=="text":
3991           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3201  def power(arg0,arg1): Line 4046  def power(arg0,arg1):
4046         """         """
4047         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4048         if testForZero(args[0]):         if testForZero(args[0]):
4049            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4050         elif testForZero(args[1]):         elif testForZero(args[1]):
4051            return numarray.ones(args[0],numarray.Float)            return numarray.ones(getShape(args[1]),numarray.Float64)
4052         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4053            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4054         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 4071  class Power_Symbol(DependendSymbol):
4071         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4072         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4073         """         """
4074         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4075         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4076         if not sh0==sh1:         if not sh0==sh1:
4077            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4078         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3244  class Power_Symbol(DependendSymbol): Line 4089  class Power_Symbol(DependendSymbol):
4089        @type format: C{str}        @type format: C{str}
4090        @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.
4091        @rtype: C{str}        @rtype: C{str}
4092        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4093        """        """
4094        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4095           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3304  def maximum(*args): Line 4149  def maximum(*args):
4149         if out==None:         if out==None:
4150            out=a            out=a
4151         else:         else:
4152            m=whereNegative(out-a)            diff=add(a,-out)
4153            out=m*a+(1.-m)*out            out=add(out,mult(wherePositive(diff),diff))
4154      return out      return out
4155        
4156  def minimum(*arg):  def minimum(*args):
4157      """      """
4158      the minimum over arguments args      the minimum over arguments args
4159    
# Line 3322  def minimum(*arg): Line 4167  def minimum(*arg):
4167         if out==None:         if out==None:
4168            out=a            out=a
4169         else:         else:
4170            m=whereNegative(out-a)            diff=add(a,-out)
4171            out=m*out+(1.-m)*a            out=add(out,mult(whereNegative(diff),diff))
4172      return out      return out
4173    
4174    def clip(arg,minval=None,maxval=None):
4175        """
4176        cuts the values of arg between minval and maxval
4177    
4178        @param arg: argument
4179        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4180        @param minval: lower range. If None no lower range is applied
4181        @type minval: C{float} or C{None}
4182        @param maxval: upper range. If None no upper range is applied
4183        @type maxval: C{float} or C{None}
4184        @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4185                 less then maxval are unchanged.
4186        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4187        @raise ValueError: if minval>maxval
4188        """
4189        if not minval==None and not maxval==None:
4190           if minval>maxval:
4191              raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4192        if minval == None:
4193            tmp=arg
4194        else:
4195            tmp=maximum(minval,arg)
4196        if maxval == None:
4197            return tmp
4198        else:
4199            return minimum(tmp,maxval)
4200    
4201        
4202  def inner(arg0,arg1):  def inner(arg0,arg1):
4203      """      """
# Line 3340  def inner(arg0,arg1): Line 4213  def inner(arg0,arg1):
4213      @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}
4214      @param arg1: second argument      @param arg1: second argument
4215      @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}
4216      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4217      @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
4218      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4219      """      """
4220      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4221      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4222      if not sh0==sh1:      if not sh0==sh1:
4223          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4224      return generalTensorProduct(arg0,arg1,offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4225    
4226    def outer(arg0,arg1):
4227        """
4228        the outer product of the two argument:
4229    
4230        out[t,s]=arg0[t]*arg1[s]
4231    
4232        where
4233    
4234            - s runs through arg0.Shape
4235            - t runs through arg1.Shape
4236    
4237        @param arg0: first argument
4238        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4239        @param arg1: second argument
4240        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4241        @return: the outer product of arg0 and arg1 at each data point
4242        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4243        """
4244        return generalTensorProduct(arg0,arg1,axis_offset=0)
4245    
4246  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4247      """      """
4248        see L{matrix_mult}
4249        """
4250        return matrix_mult(arg0,arg1)
4251    
4252    def matrix_mult(arg0,arg1):
4253        """
4254      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4255    
4256      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4257    
4258            or      or
4259    
4260      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4261    
4262      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.
4263    
4264      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4265      @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 4269  def matrixmult(arg0,arg1):
4269      @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
4270      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4271      """      """
4272      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4273      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4274      if not len(sh0)==2 :      if not len(sh0)==2 :
4275          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4276      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4277          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4278      return generalTensorProduct(arg0,arg1,offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4279    
4280  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4281      """      """
4282      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  
4283      """      """
4284      return generalTensorProduct(arg0,arg1,offset=0)      return tensor_mult(arg0,arg1)
4285    
4286    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4287      """      """
4288      the tensor product of the two argument:      the tensor product of the two argument:
   
4289            
4290      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4291    
4292      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4293    
4294                   or      or
4295    
4296      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4297    
# Line 3415  def tensormult(arg0,arg1): Line 4300  def tensormult(arg0,arg1):
4300    
4301      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]
4302                                
4303                   or      or
4304    
4305      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]
4306    
4307                   or      or
4308    
4309      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]
4310    
4311      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  
4312      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.
4313    
4314      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4315      @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 4318  def tensormult(arg0,arg1):
4318      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4319      @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
4320      """      """
4321      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4322      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4323      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4324         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4325      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):
4326         return generalTensorProduct(arg0,arg1,offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4327      else:      else:
4328          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4329    
4330  def generalTensorProduct(arg0,arg1,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4331      """      """
4332      generalized tensor product      generalized tensor product
4333    
4334      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4335    
4336      where s runs through arg0.Shape[:arg0.Rank-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[offset:]  
4337    
4338      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]
4339      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4340            - t runs through arg1.Shape[axis_offset:]
4341    
4342      @param arg0: first argument      @param arg0: first argument
4343      @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}
4344      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4345      @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}
4346      @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.
4347      @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 4351  def generalTensorProduct(arg0,arg1,offse
4351      # 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
4352      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4353         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4354             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4355         else:         else:
4356             if not arg0.shape[arg0.rank-offset:]==arg1.shape[:offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4357                 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)
4358             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4359             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4360             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
4361             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4362             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4363             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4364             for i in sh1[:offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4365             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4366             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4367             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4368             for i0 in range(d0):             for i0 in range(d0):
4369                      for i1 in range(d1):                      for i1 in range(d1):
4370                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
4371             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4372             return out             return out
4373      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4374         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4375             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4376         else:         else:
4377             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)
4378      else:            else:      
4379         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4380                                    
4381  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4382     """     """
4383     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4384     """     """
4385     def __init__(self,arg0,arg1,offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4386         """         """
4387         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4388    
4389         @param arg0: numerator         @param arg0: first argument
4390         @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}.
4391         @param arg1: denominator         @param arg1: second argument
4392         @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}.
4393         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4394         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4395         """         """
4396         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4397         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4398         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4399         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4400         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4401         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4402         if not sh01==sh10:         if not sh01==sh10:
4403             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)
4404         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])
4405    
4406     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4407        """        """
# Line 3529  class GeneralTensorProduct_Symbol(Depend Line 4413  class GeneralTensorProduct_Symbol(Depend
4413        @type format: C{str}        @type format: C{str}
4414        @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.
4415        @rtype: C{str}        @rtype: C{str}
4416        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4417        """        """
4418        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4419           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])
4420        else:        else:
4421           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)
4422    
# Line 3557  class GeneralTensorProduct_Symbol(Depend Line 4441  class GeneralTensorProduct_Symbol(Depend
4441           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4442           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4443    
4444  def escript_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4445      "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!!!"
4446      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4447      shape0=arg0.getShape()[:arg0.getRank()-offset]  
4448      shape01=arg0.getShape()[arg0.getRank()-offset:]  def transposed_matrix_mult(arg0,arg1):
4449      shape10=arg1.getShape()[:offset]      """
4450      shape1=arg1.getShape()[offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4451      if not shape01==shape10:  
4452          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]
4453    
4454      # create return value:      or
4455      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())  
4456      #      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4457      s0=[[]]  
4458      for k in shape0:      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4459            s=[]  
4460            for j in s0:      The first dimension of arg0 and arg1 must match.
4461                  for i in range(k): s.append(j+[slice(i,i)])  
4462            s0=s      @param arg0: first argument of rank 2
4463      s1=[[]]      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4464      for k in shape1:      @param arg1: second argument of at least rank 1
4465            s=[]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4466            for j in s1:      @return: the product of the transposed of arg0 and arg1 at each data point
4467                  for i in range(k): s.append(j+[slice(i,i)])      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4468            s1=s      @raise ValueError: if the shapes of the arguments are not appropriate
4469      s01=[[]]      """
4470      for k in shape01:      sh0=getShape(arg0)
4471            s=[]      sh1=getShape(arg1)
4472            for j in s01:      if not len(sh0)==2 :
4473                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"first argument must have rank 2"
4474            s01=s      if not len(sh1)==2 and not len(sh1)==1:
4475            raise ValueError,"second argument must have rank 1 or 2"
4476      for i0 in s0:      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4477         for i1 in s1:  
4478           s=escript.Scalar(0.,arg0.getFunctionSpace())  def transposed_tensor_mult(arg0,arg1):
4479           for i01 in s01:      """
4480              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      the tensor product of the transposed of the first and the second argument
4481           out.__setitem__(tuple(i0+i1),s)      
4482      return out      for arg0 of rank 2 this is
4483    
4484        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4485    
4486        or
4487    
4488        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4489    
4490      
4491        and for arg0 of rank 4 this is
4492    
4493        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4494                  
4495        or
4496    
4497        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4498    
4499        or
4500    
4501        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4502    
4503        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4504        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4505    
4506        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4507    
4508        @param arg0: first argument of rank 2 or 4
4509        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4510        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4511        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4512        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4513        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4514        """
4515        sh0=getShape(arg0)
4516        sh1=getShape(arg1)
4517        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4518           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4519        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4520           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4521        else:
4522            raise ValueError,"first argument must have rank 2 or 4"
4523    
4524    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4525        """
4526        generalized tensor product of transposed of arg0 and arg1:
4527    
4528        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4529    
4530        where
4531    
4532            - s runs through arg0.Shape[axis_offset:]
4533            - r runs trough arg0.Shape[:axis_offset]
4534            - t runs through arg1.Shape[axis_offset:]
4535    
4536        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4537        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4538    
4539        @param arg0: first argument
4540        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4541        @param arg1: second argument
4542        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4543        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4544        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4545        """
4546        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4547        arg0,arg1=matchType(arg0,arg1)
4548        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4549        if isinstance(arg0,numarray.NumArray):
4550           if isinstance(arg1,Symbol):
4551               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4552           else:
4553               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4554                   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)
4555               arg0_c=arg0.copy()
4556               arg1_c=arg1.copy()
4557               sh0,sh1=arg0.shape,arg1.shape
4558               d0,d1,d01=1,1,1
4559               for i in sh0[axis_offset:]: d0*=i
4560               for i in sh1[axis_offset:]: d1*=i
4561               for i in sh0[:axis_offset]: d01*=i
4562               arg0_c.resize((d01,d0))
4563               arg1_c.resize((d01,d1))
4564               out=numarray.zeros((d0,d1),numarray.Float64)
4565               for i0 in range(d0):
4566                        for i1 in range(d1):
4567                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4568               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4569               return out
4570        elif isinstance(arg0,escript.Data):
4571           if isinstance(arg1,Symbol):
4572               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4573           else:
4574               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4575        else:      
4576           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4577                    
4578    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4579       """
4580       Symbol representing the general tensor product of the transposed of arg0 and arg1
4581       """
4582       def __init__(self,arg0,arg1,axis_offset=0):
4583           """
4584           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4585    
4586           @param arg0: first argument
4587           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4588           @param arg1: second argument
4589           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4590           @raise ValueError: inconsistent dimensions of arguments.
4591           @note: if both arguments have a spatial dimension, they must equal.
4592           """
4593           sh_arg0=getShape(arg0)
4594           sh_arg1=getShape(arg1)
4595           sh01=sh_arg0[:axis_offset]
4596           sh10=sh_arg1[:axis_offset]
4597           sh0=sh_arg0[axis_offset:]
4598           sh1=sh_arg1[axis_offset:]
4599           if not sh01==sh10:
4600               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)
4601           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4602    
4603       def getMyCode(self,argstrs,format="escript"):
4604          """
4605          returns a program code that can be used to evaluate the symbol.
4606    
4607          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4608          @type argstrs: C{list} of length 2 of C{str}.
4609          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4610          @type format: C{str}
4611          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4612          @rtype: C{str}
4613          @raise NotImplementedError: if the requested format is not available
4614          """
4615          if format=="escript" or format=="str" or format=="text":
4616             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4617          else:
4618             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4619    
4620       def substitute(self,argvals):
4621          """
4622          assigns new values to symbols in the definition of the symbol.
4623          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4624    
4625          @param argvals: new values assigned to symbols
4626          @type argvals: C{dict} with keywords of type L{Symbol}.
4627          @return: result of the substitution process. Operations are executed as much as possible.
4628          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4629          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4630          """
4631          if argvals.has_key(self):
4632             arg=argvals[self]
4633             if self.isAppropriateValue(arg):
4634                return arg
4635             else:
4636                raise TypeError,"%s: new value is not appropriate."%str(self)
4637          else:
4638             args=self.getSubstitutedArguments(argvals)
4639             return generalTransposedTensorProduct(args[0],args[1],args[2])
4640    
4641    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4642        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4643        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4644    
4645    def matrix_transposed_mult(arg0,arg1):
4646        """
4647        matrix-transposed(matrix) product of the two argument:
4648    
4649        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4650    
4651        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4652    
4653        The last dimensions of arg0 and arg1 must match.
4654    
4655        @param arg0: first argument of rank 2
4656        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4657        @param arg1: second argument of rank 2
4658        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4659        @return: the product of arg0 and the transposed of arg1 at each data point
4660        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4661        @raise ValueError: if the shapes of the arguments are not appropriate
4662        """
4663        sh0=getShape(arg0)
4664        sh1=getShape(arg1)
4665        if not len(sh0)==2 :
4666            raise ValueError,"first argument must have rank 2"
4667        if not len(sh1)==2 and not len(sh1)==1:
4668            raise ValueError,"second argument must have rank 1 or 2"
4669        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4670    
4671    def tensor_transposed_mult(arg0,arg1):
4672        """
4673        the tensor product of the first and the transpose of the second argument
4674        
4675        for arg0 of rank 2 this is
4676    
4677        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4678    
4679        and for arg0 of rank 4 this is
4680    
4681        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4682                  
4683        or
4684    
4685        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4686    
4687        In the first case the the second dimension of arg0 and arg1 must match and  
4688        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4689    
4690        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4691    
4692        @param arg0: first argument of rank 2 or 4
4693        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4694        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4695        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4696        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4697        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4698        """
4699        sh0=getShape(arg0)
4700        sh1=getShape(arg1)
4701        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4702           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4703        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4704           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4705        else:
4706            raise ValueError,"first argument must have rank 2 or 4"
4707    
4708    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4709        """
4710        generalized tensor product of transposed of arg0 and arg1:
4711    
4712        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4713    
4714        where
4715    
4716            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4717            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4718            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4719    
4720        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4721        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4722    
4723        @param arg0: first argument
4724        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4725        @param arg1: second argument
4726        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4727        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4728        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4729        """
4730        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4731        arg0,arg1=matchType(arg0,arg1)
4732        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4733        if isinstance(arg0,numarray.NumArray):
4734           if isinstance(arg1,Symbol):
4735               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4736           else:
4737               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4738                   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)
4739               arg0_c=arg0.copy()
4740               arg1_c=arg1.copy()
4741               sh0,sh1=arg0.shape,arg1.shape
4742               d0,d1,d01=1,1,1
4743               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4744               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4745               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4746               arg0_c.resize((d0,d01))
4747               arg1_c.resize((d1,d01))
4748               out=numarray.zeros((d0,d1),numarray.Float64)
4749               for i0 in range(d0):
4750                        for i1 in range(d1):
4751                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4752               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4753               return out
4754        elif isinstance(arg0,escript.Data):
4755           if isinstance(arg1,Symbol):
4756               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4757           else:
4758               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4759        else:      
4760           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4761                    
4762    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4763       """
4764       Symbol representing the general tensor product of arg0 and the transpose of arg1
4765       """
4766       def __init__(self,arg0,arg1,axis_offset=0):
4767           """
4768           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4769    
4770           @param arg0: first argument
4771           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4772           @param arg1: second argument
4773           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4774           @raise ValueError: inconsistent dimensions of arguments.
4775           @note: if both arguments have a spatial dimension, they must equal.
4776           """
4777           sh_arg0=getShape(arg0)
4778           sh_arg1=getShape(arg1)
4779           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4780           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4781           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4782           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4783           if not sh01==sh10:
4784               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)
4785           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4786    
4787       def getMyCode(self,argstrs,format="escript"):
4788          """
4789          returns a program code that can be used to evaluate the symbol.
4790    
4791          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4792          @type argstrs: C{list} of length 2 of C{str}.
4793          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4794          @type format: C{str}
4795          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4796          @rtype: C{str}
4797          @raise NotImplementedError: if the requested format is not available
4798          """
4799          if format=="escript" or format=="str" or format=="text":
4800             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4801          else:
4802             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4803    
4804       def substitute(self,argvals):
4805          """
4806          assigns new values to symbols in the definition of the symbol.
4807          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4808    
4809          @param argvals: new values assigned to symbols
4810          @type argvals: C{dict} with keywords of type L{Symbol}.
4811          @return: result of the substitution process. Operations are executed as much as possible.
4812          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4813          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4814          """
4815          if argvals.has_key(self):
4816             arg=argvals[self]
4817             if self.isAppropriateValue(arg):
4818                return arg
4819             else:
4820                raise TypeError,"%s: new value is not appropriate."%str(self)
4821          else:
4822             args=self.getSubstitutedArguments(argvals)
4823             return generalTensorTransposedProduct(args[0],args[1],args[2])
4824    
4825    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4826        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4827        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4828    
4829  #=========================================================  #=========================================================
4830  #   some little helpers  #  functions dealing with spatial dependency
4831  #=========================================================  #=========================================================
4832  def grad(arg,where=None):  def grad(arg,where=None):
4833      """      """
4834      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4835    
4836        If C{g} is the returned object, then
4837    
4838          - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4839          - 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.
4840          - 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.
4841          - 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.
4842    
4843      @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.
4844                    to be calculated.      @type arg: L{escript.Data} or L{Symbol}
4845      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
4846                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4847        @type where: C{None} or L{escript.FunctionSpace}
4848        @return: gradient of arg.
4849        @rtype: L{escript.Data} or L{Symbol}
4850      """      """
4851      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4852         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3617  def grad(arg,where=None): Line 4856  def grad(arg,where=None):
4856         else:         else:
4857            return arg._grad(where)            return arg._grad(where)
4858      else:      else:
4859        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4860    
4861    class Grad_Symbol(DependendSymbol):
4862       """
4863       L{Symbol} representing the result of the gradient operator
4864       """
4865       def __init__(self,arg,where=None):
4866          """
4867          initialization of gradient L{Symbol} with argument arg
4868          @param arg: argument of function
4869          @type arg: L{Symbol}.
4870          @param where: FunctionSpace in which the gradient will be calculated.
4871                      If not present or C{None} an appropriate default is used.
4872          @type where: C{None} or L{escript.FunctionSpace}
4873          """
4874          d=arg.getDim()
4875          if d==None:
4876             raise ValueError,"argument must have a spatial dimension"
4877          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4878    
4879       def getMyCode(self,argstrs,format="escript"):
4880          """
4881          returns a program code that can be used to evaluate the symbol.
4882    
4883          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4884          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4885          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4886          @type format: C{str}
4887          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4888          @rtype: C{str}
4889          @raise NotImplementedError: if the requested format is not available
4890          """
4891          if format=="escript" or format=="str"  or format=="text":
4892             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4893          else:
4894             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4895    
4896       def substitute(self,argvals):
4897          """
4898          assigns new values to symbols in the definition of the symbol.
4899          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4900    
4901          @param argvals: new values assigned to symbols
4902          @type argvals: C{dict} with keywords of type L{Symbol}.
4903          @return: result of the substitution process. Operations are executed as much as possible.
4904          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4905          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4906          """
4907          if argvals.has_key(self):
4908             arg=argvals[self]
4909             if self.isAppropriateValue(arg):
4910                return arg
4911             else:
4912                raise TypeError,"%s: new value is not appropriate."%str(self)
4913          else:
4914             arg=self.getSubstitutedArguments(argvals)
4915             return grad(arg[0],where=arg[1])
4916    
4917       def diff(self,arg):
4918          """
4919          differential of this object
4920    
4921          @param arg: the derivative is calculated with respect to arg
4922          @type arg: L{escript.Symbol}
4923          @return: derivative with respect to C{arg}
4924          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4925          """
4926          if arg==self:
4927             return identity(self.getShape())
4928          else:
4929             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4930    
4931  def integrate(arg,where=None):  def integrate(arg,where=None):
4932      """      """
4933      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}
4934      its domain.      before integration.
4935    
4936      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4937        @type arg: L{escript.Data} or L{Symbol}
4938      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
4939                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4940        @type where: C{None} or L{escript.FunctionSpace}
4941        @return: integral of arg.
4942        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4943      """      """
4944      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):  
4945         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4946      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4947         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 4950  def integrate(arg,where=None):
4950         else:         else:
4951            return arg._integrate()            return arg._integrate()
4952      else:      else:
4953        raise TypeError,"integrate: Unknown argument type."         arg2=escript.Data(arg,where)
4954           if arg2.getRank()==0:
4955              return arg2._integrate()[0]
4956           else:
4957              return arg2._integrate()
4958    
4959    class Integrate_Symbol(DependendSymbol):
4960       """
4961       L{Symbol} representing the result of the spatial integration operator
4962       """
4963       def __init__(self,arg,where=None):
4964          """
4965          initialization of integration L{Symbol} with argument arg
4966          @param arg: argument of the integration
4967          @type arg: L{Symbol}.
4968          @param where: FunctionSpace in which the integration will be calculated.
4969                      If not present or C{None} an appropriate default is used.
4970          @type where: C{None} or L{escript.FunctionSpace}
4971          """
4972          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4973    
4974       def getMyCode(self,argstrs,format="escript"):
4975          """
4976          returns a program code that can be used to evaluate the symbol.
4977    
4978          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4979          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4980          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4981          @type format: C{str}
4982          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4983          @rtype: C{str}
4984          @raise NotImplementedError: if the requested format is not available
4985          """
4986          if format=="escript" or format=="str"  or format=="text":
4987             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4988          else:
4989             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4990    
4991       def substitute(self,argvals):
4992          """
4993          assigns new values to symbols in the definition of the symbol.
4994          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4995    
4996          @param argvals: new values assigned to symbols
4997          @type argvals: C{dict} with keywords of type L{Symbol}.
4998          @return: result of the substitution process. Operations are executed as much as possible.
4999          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
5000          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
5001          """
5002          if argvals.has_key(self):
5003             arg=argvals[self]
5004             if self.isAppropriateValue(arg):
5005                return arg
5006             else:
5007                raise TypeError,"%s: new value is not appropriate."%str(self)
5008          else:
5009             arg=self.getSubstitutedArguments(argvals)
5010             return integrate(arg[0],where=arg[1])
5011    
5012       def diff(self,arg):
5013          """
5014          differential of this object
5015    
5016          @param arg: the derivative is calculated with respect to arg
5017          @type arg: L{escript.Symbol}
5018          @return: derivative with respect to C{arg}
5019          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
5020          """
5021          if arg==self:
5022             return identity(self.getShape())
5023          else:
5024             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5025    
5026    
5027  def interpolate(arg,where):  def interpolate(arg,where):
5028      """      """
5029      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where. If the argument C{arg} has the requested function space
5030        C{where} no interpolation is performed and C{arg} is returned.
5031    
5032      @param arg:    interpolant      @param arg: interpolant
5033      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
5034        @param where: FunctionSpace to be interpolated to
5035        @type where: L{escript.FunctionSpace}
5036        @return: interpolated argument
5037        @rtype: C{escript.Data} or L{Symbol}
5038      """      """
5039      if testForZero(arg):      if isinstance(arg,Symbol):
5040        return 0         return Interpolate_Symbol(arg,where)
5041      elif isinstance(arg,Symbol):      elif isinstance(arg,escript.Data):
5042         return Interpolated_Symbol(arg,where)         if arg.isEmpty():
5043              return arg
5044           elif where == arg.getFunctionSpace():
5045              return arg
5046           else:
5047              return escript.Data(arg,where)
5048      else:      else:
5049         return escript.Data(arg,where)         return escript.Data(arg,where)
5050    
5051    class Interpolate_Symbol(DependendSymbol):
5052       """
5053       L{Symbol} representing the result of the interpolation operator
5054       """
5055       def __init__(self,arg,where):
5056          """
5057          initialization of interpolation L{Symbol} with argument arg
5058          @param arg: argument of the interpolation
5059          @type arg: L{Symbol}.
5060          @param where: FunctionSpace into which the argument is interpolated.
5061          @type where: L{escript.FunctionSpace}
5062          """
5063          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
5064    
5065       def getMyCode(self,argstrs,format="escript"):
5066          """
5067          returns a program code that can be used to evaluate the symbol.
5068    
5069          @param argstrs: gives for each argument a string representing the argument for the evaluation.
5070          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
5071          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
5072          @type format: C{str}
5073          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
5074          @rtype: C{str}
5075          @raise NotImplementedError: if the requested format is not available
5076          """
5077          if format=="escript" or format=="str"  or format=="text":
5078             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
5079          else:
5080             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
5081    
5082       def substitute(self,argvals):
5083          """
5084          assigns new values to symbols in the definition of the symbol.
5085          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
5086    
5087          @param argvals: new values assigned to symbols
5088          @type argvals: C{dict} with keywords of type L{Symbol}.
5089          @return: result of the substitution process. Operations are executed as much as possible.
5090          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
5091          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
5092          """
5093          if argvals.has_key(self):
5094             arg=argvals[self]
5095             if self.isAppropriateValue(arg):
5096                return arg
5097             else:
5098                raise TypeError,"%s: new value is not appropriate."%str(self)
5099          else:
5100             arg=self.getSubstitutedArguments(argvals)
5101             return interpolate(arg[0],where=arg[1])
5102    
5103       def diff(self,arg):
5104          """
5105          differential of this object
5106    
5107          @param arg: the derivative is calculated with respect to arg
5108          @type arg: L{escript.Symbol}
5109          @return: derivative with respect to C{arg}
5110          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
5111          """
5112          if arg==self:
5113             return identity(self.getShape())
5114          else:
5115             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
5116    
5117    
5118  def div(arg,where=None):  def div(arg,where=None):
5119      """      """
5120      Returns the divergence of arg at where.      returns the divergence of arg at where.
5121    
5122      @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.
5123                    be calculated.      @type arg: L{escript.Data} or L{Symbol}
5124      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the divergence will be calculated.
5125                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
5126        @type where: C{None} or L{escript.FunctionSpace}
5127        @return: divergence of arg.
5128        @rtype: L{escript.Data} or L{Symbol}
5129      """      """
5130      g=grad(arg,where)      if isinstance(arg,Symbol):
5131      return trace(g,axis0=g.getRank()-2,axis1=g.getRank()-1)          dim=arg.getDim()
5132        elif isinstance(arg,escript.Data):
5133            dim=arg.getDomain().getDim()
5134        else:
5135            raise TypeError,"div: argument type not supported"
5136        if not arg.getShape()==(dim,):
5137          raise ValueError,"div: expected shape is (%s,)"%dim
5138        return trace(grad(arg,where))
5139    
5140  def jump(arg):  def jump(arg,domain=None):
5141      """      """
5142      Returns the jump of arg across a continuity.      returns the jump of arg across the continuity of the domain
5143    
5144      @param arg:   Data object representing the function which gradient      @param arg: argument
5145                    to be calculated.      @type arg: L{escript.Data} or L{Symbol}
5146        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5147                       the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5148        @type domain: C{None} or L{escript.Domain}
5149        @return: jump of arg
5150        @rtype: L{escript.Data} or L{Symbol}
5151      """      """
5152      d=arg.getDomain()      if domain==None: domain=arg.getDomain()
5153      return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
5154    
5155  #=============================  def L2(arg):
5156  #      """
5157  # wrapper for various functions: if the argument has attribute the function name      returns the L2 norm of arg at where
5158  # as an argument it calls the corresponding methods. Otherwise the corresponding      
5159  # numarray function is called.      @param arg: function which L2 to be calculated.
5160        @type arg: L{escript.Data} or L{Symbol}
5161  # functions involving the underlying Domain:      @return: L2 norm of arg.
5162        @rtype: L{float} or L{Symbol}
5163        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5164        """
5165        return sqrt(integrate(inner(arg,arg)))
5166    
5167  def transpose(arg,axis=None):  def getClosestValue(arg,origin=0):
5168        """
5169        returns the value in arg which is closest to origin
5170        
5171        @param arg: function
5172        @type arg: L{escript.Data} or L{Symbol}
5173        @param origin: reference value
5174        @type origin: C{float} or L{escript.Data}
5175        @return: value in arg closest to origin
5176        @rtype: L{numarray.NumArray} or L{Symbol}
5177      """      """
5178      Returns the transpose of the Data object arg.      return arg.getValueOfGlobalDataPoint(*(length(arg-origin).minGlobalDataPoint()))
5179    
5180      @param arg:  def normalize(arg,zerolength=0):
5181      """      """
5182      if axis==None:      returns normalized version of arg (=arg/length(arg))
5183         r=0      
5184         if hasattr(arg,"getRank"): r=arg.getRank()      @param arg: function
5185         if hasattr(arg,"rank"): r=arg.rank      @type arg: L{escript.Data} or L{Symbol}
5186         axis=r/2      @param zerolength: realitive tolerance for arg == 0.
5187      if isinstance(arg,Symbol):      @type zerolength: C{float}
5188         return Transpose_Symbol(arg,axis=r)      @return: normalized arg where arg is non zero and zero elsewhere
5189      if isinstance(arg,escript.Data):      @rtype: L{escript.Data} or L{Symbol}
5190         # hack for transpose      """
5191         r=arg.getRank()