/[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 429 by gross, Wed Jan 11 05:53:40 2006 UTC revision 2142 by jfenwick, Tue Dec 9 06:22:57 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$"  
 __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
45  # def pokeShape(arg):  from esys.escript import listEscriptParams
 # def pokeDim(arg):  
 # def commonShape(arg0,arg1):  
 # def commonDim(*args):  
 # def testForZero(arg):  
 # def matchType(arg0=0.,arg1=0.):  
 # def matchShape(arg0,arg1):  
   
 # def transpose(arg,axis=None):  
 # def trace(arg,axis0=0,axis1=1):  
 # def reorderComponents(arg,index):  
   
 # def integrate(arg,where=None):  
 # def interpolate(arg,where):  
 # def div(arg,where=None):  
 # def grad(arg,where=None):  
   
 #  
 # slicing: get  
 #          set  
 #  
 # and derivatives  
46    
47  #=========================================================  #=========================================================
48  #   some helpers:  #   some helpers:
49  #=========================================================  #=========================================================
50    def getEpsilon():
51         return escript.getMachinePrecision()
52    EPSILON=getEpsilon()
53    
54    def getMaxFloat():
55         return escript.getMaxFloat()
56    DBLE_MAX=getMaxFloat()
57    
58    
59    def getTagNames(domain):
60        """
61        returns a list of the tag names used by the domain
62    
63        
64        @param domain: a domain object
65        @type domain: L{escript.Domain}
66        @return: a list of the tag name used by the domain.
67        @rtype: C{list} of C{str}
68        """
69        return [n.strip() for n in domain.showTagNames().split(",") ]
70    
71    def insertTagNames(domain,**kwargs):
72        """
73        inserts tag names into the domain
74    
75        @param domain: a domain object
76        @type domain: C{escript.Domain}
77        @keyword <tag_name>: tag key assigned to <tag_name>
78        @type <tag_name>: C{int}
79        """
80        for  k in kwargs:
81             domain.setTagMap(k,kwargs[k])
82    
83    def insertTaggedValues(target,**kwargs):
84        """
85        inserts tagged values into the tagged using tag names
86    
87        @param target: data to be filled by tagged values
88        @type target: L{escript.Data}
89        @keyword <tag_name>: value to be used for <tag_name>
90        @type <tag_name>: C{float} or {numarray.NumArray}
91        @return: C{target}
92        @rtype: L{escript.Data}
93        """
94        for k in kwargs:
95            target.setTaggedValue(k,kwargs[k])
96        return target
97    
98  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
99      """      """
100      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.
101    
102      Example:      Example::
103    
104         tmp=Scalar(..)         tmp=Scalar(..)
105         v=Vector(..)         v=Vector(..)
106         saveVTK("solution.xml",temperature=tmp,velovity=v)         saveVTK("solution.xml",temperature=tmp,velocity=v)
107    
108      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"
109    
110      @param filename: file name of the output file      @param filename: file name of the output file
111      @type filename: C{str}      @type filename: C{str}
# Line 81  def saveVTK(filename,domain=None,**data) Line 115  def saveVTK(filename,domain=None,**data)
115      @type <name>: L{Data} object.      @type <name>: L{Data} object.
116      @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.
117      """      """
118      if domain==None:      new_data={}
119         for i in data.keys():      for n,d in data.items():
120            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
121                fs=d.getFunctionSpace()
122                domain2=fs.getDomain()
123                if fs == escript.Solution(domain2):
124                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
125                elif fs == escript.ReducedSolution(domain2):
126                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
127                else:
128                   new_data[n]=d
129                if domain==None: domain=domain2
130      if domain==None:      if domain==None:
131          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
132      else:      domain.saveVTK(filename,new_data)
         domain.saveVTK(filename,data)  
133    
134  def saveDX(filename,domain=None,**data):  def saveDX(filename,domain=None,**data):
135      """      """
136      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.
137    
138      Example:      Example::
139    
140         tmp=Scalar(..)         tmp=Scalar(..)
141         v=Vector(..)         v=Vector(..)
142         saveDX("solution.dx",temperature=tmp,velovity=v)         saveDX("solution.dx",temperature=tmp,velocity=v)
143    
144      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".
145    
146      @param filename: file name of the output file      @param filename: file name of the output file
147      @type filename: C{str}      @type filename: C{str}
# Line 109  def saveDX(filename,domain=None,**data): Line 151  def saveDX(filename,domain=None,**data):
151      @type <name>: L{Data} object.      @type <name>: L{Data} object.
152      @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.
153      """      """
154      if domain==None:      new_data={}
155         for i in data.keys():      for n,d in data.items():
156            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
157                fs=d.getFunctionSpace()
158                domain2=fs.getDomain()
159                if fs == escript.Solution(domain2):
160                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
161                elif fs == escript.ReducedSolution(domain2):
162                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
163                elif fs == escript.ContinuousFunction(domain2):
164                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
165                else:
166                   new_data[n]=d
167                if domain==None: domain=domain2
168      if domain==None:      if domain==None:
169          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
170      else:      domain.saveDX(filename,new_data)
         domain.saveDX(filename,data)  
171    
172  def kronecker(d=3):  def kronecker(d=3):
173     """     """
174     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
175    
176     @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
177     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
178     @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
179     @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}  
180     """     """
181     return identityTensor(d)     return identityTensor(d)
182    
# Line 140  def identity(shape=()): Line 191  def identity(shape=()):
191     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
192     """     """
193     if len(shape)>0:     if len(shape)>0:
194        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
195        if len(shape)==1:        if len(shape)==1:
196            for i0 in range(shape[0]):            for i0 in range(shape[0]):
197               out[i0,i0]=1.               out[i0,i0]=1.
   
198        elif len(shape)==2:        elif len(shape)==2:
199            for i0 in range(shape[0]):            for i0 in range(shape[0]):
200               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 160  def identityTensor(d=3): Line 210  def identityTensor(d=3):
210     return the dxd identity matrix     return the dxd identity matrix
211    
212     @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
213     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
214     @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
215     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
216     """     """
217     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
218        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
219     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
220           return identity((d.getDim(),))
221       else:
222           return identity((d,))
223    
224  def identityTensor4(d=3):  def identityTensor4(d=3):
225     """     """
# Line 175  def identityTensor4(d=3): Line 228  def identityTensor4(d=3):
228     @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
229     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
230     @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
231     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
232     """     """
233     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
234        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
235     return identity((d,d))     elif isinstance(d,escript.Domain):
236           return identity((d.getDim(),d.getDim()))
237       else:
238           return identity((d,d))
239    
240  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
241     """     """
# Line 188  def unitVector(i=0,d=3): Line 244  def unitVector(i=0,d=3):
244     @param i: index     @param i: index
245     @type i: C{int}     @type i: C{int}
246     @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
247     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
248     @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
249     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
250     """     """
251     return kronecker(d)[i]     return kronecker(d)[i]
252    
# Line 246  def inf(arg): Line 302  def inf(arg):
302    
303      @param arg: argument      @param arg: argument
304      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
305      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
306      @rtype: C{float}      @rtype: C{float}
307      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
308      """      """
# Line 265  def inf(arg): Line 321  def inf(arg):
321  #=========================================================================  #=========================================================================
322  #   some little helpers  #   some little helpers
323  #=========================================================================  #=========================================================================
324  def pokeShape(arg):  def getRank(arg):
325        """
326        identifies the rank of its argument
327    
328        @param arg: a given object
329        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
330        @return: the rank of the argument
331        @rtype: C{int}
332        @raise TypeError: if type of arg cannot be processed
333        """
334    
335        if isinstance(arg,numarray.NumArray):
336            return arg.rank
337        elif isinstance(arg,escript.Data):
338            return arg.getRank()
339        elif isinstance(arg,float):
340            return 0
341        elif isinstance(arg,int):
342            return 0
343        elif isinstance(arg,Symbol):
344            return arg.getRank()
345        else:
346          raise TypeError,"getShape: cannot identify shape"
347    def getShape(arg):
348      """      """
349      identifies the shape of its argument      identifies the shape of its argument
350    
# Line 287  def pokeShape(arg): Line 366  def pokeShape(arg):
366      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
367          return arg.getShape()          return arg.getShape()
368      else:      else:
369        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
370    
371  def pokeDim(arg):  def pokeDim(arg):
372      """      """
# Line 310  def commonShape(arg0,arg1): Line 389  def commonShape(arg0,arg1):
389      """      """
390      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.
391    
392      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
393      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
394      @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.
395      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
396      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
397      """      """
398      sh0=pokeShape(arg0)      sh0=getShape(arg0)
399      sh1=pokeShape(arg1)      sh1=getShape(arg1)
400      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
401         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
402               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
# Line 335  def commonDim(*args): Line 414  def commonDim(*args):
414      """      """
415      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.
416    
417      @param *args: given objects      @param args: given objects
418      @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
419               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
420      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 357  def testForZero(arg): Line 436  def testForZero(arg):
436    
437      @param arg: a given object      @param arg: a given object
438      @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}
439      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
440      @rtype : C{bool}      @rtype: C{bool}
441      """      """
442      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
443         return not Lsup(arg)>0.         return not Lsup(arg)>0.
# Line 391  def matchType(arg0=0.,arg1=0.): Line 470  def matchType(arg0=0.,arg1=0.):
470         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
471            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
472         elif isinstance(arg1,float):         elif isinstance(arg1,float):
473            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
474         elif isinstance(arg1,int):         elif isinstance(arg1,int):
475            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
476         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
477            pass            pass
478         else:         else:
# Line 417  def matchType(arg0=0.,arg1=0.): Line 496  def matchType(arg0=0.,arg1=0.):
496         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
497            pass            pass
498         elif isinstance(arg1,float):         elif isinstance(arg1,float):
499            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
500         elif isinstance(arg1,int):         elif isinstance(arg1,int):
501            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
502         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
503            pass            pass
504         else:         else:
505            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
506      elif isinstance(arg0,float):      elif isinstance(arg0,float):
507         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
508            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
509         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
510            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
511         elif isinstance(arg1,float):         elif isinstance(arg1,float):
512            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
513            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
514         elif isinstance(arg1,int):         elif isinstance(arg1,int):
515            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
516            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
517         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
518            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
519         else:         else:
520            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
521      elif isinstance(arg0,int):      elif isinstance(arg0,int):
522         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
523            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
524         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
525            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
526         elif isinstance(arg1,float):         elif isinstance(arg1,float):
527            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
528            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
529         elif isinstance(arg1,int):         elif isinstance(arg1,int):
530            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
531            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
532         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
533            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
534         else:         else:
535            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
536      else:      else:
# Line 461  def matchType(arg0=0.,arg1=0.): Line 540  def matchType(arg0=0.,arg1=0.):
540    
541  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
542      """      """
543            return representations of arg0 amd arg1 which ahve the same shape
   
     If shape is not given the shape "largest" shape of args is used.  
544    
545      @param args: a given ob      @param arg0: a given object
546      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
547      @return: True if the argument is identical to zero.      @param arg1: a given object
548      @rtype: C{list} of C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
549        @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
550        @rtype: C{tuple}
551      """      """
552      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
553      sh0=pokeShape(arg0)      sh0=getShape(arg0)
554      sh1=pokeShape(arg1)      sh1=getShape(arg1)
555      if len(sh0)<len(sh):      if len(sh0)<len(sh):
556         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
557      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
558         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
559      else:      else:
560         return arg0,arg1         return arg0,arg1
561  #=========================================================  #=========================================================
# Line 496  class Symbol(object): Line 575  class Symbol(object):
575         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
576         symbols or any other object.         symbols or any other object.
577    
578         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
579         @type arg: C{list}         @type args: C{list}
580         @param shape: the shape         @param shape: the shape
581         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
582         @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
# Line 540  class Symbol(object): Line 619  class Symbol(object):
619         """         """
620         the shape of the symbol.         the shape of the symbol.
621    
622         @return : the shape of the symbol.         @return: the shape of the symbol.
623         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
624         """         """
625         return self.__shape         return self.__shape
# Line 549  class Symbol(object): Line 628  class Symbol(object):
628         """         """
629         the spatial dimension         the spatial dimension
630    
631         @return : the spatial dimension         @return: the spatial dimension
632         @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.
633         """         """
634         return self.__dim         return self.__dim
# Line 573  class Symbol(object): Line 652  class Symbol(object):
652         """         """
653         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.
654    
655         @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].
656         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
657         @rtype: C{list} of objects         @rtype: C{list} of objects
658         @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.
659         """         """
660         out=[]         out=[]
661         for a in self.getArgument():         for a in self.getArgument():
# Line 600  class Symbol(object): Line 679  class Symbol(object):
679            if isinstance(a,Symbol):            if isinstance(a,Symbol):
680               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
681            else:            else:
682                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
683                if len(s)>0:                if len(s)>0:
684                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
685                else:                else:
686                   out.append(a)                   out.append(a)
687         return out         return out
# Line 692  class Symbol(object): Line 771  class Symbol(object):
771         else:         else:
772            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
773            if len(s)>0:            if len(s)>0:
774               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
775            else:            else:
776               return 0.               return 0.
777    
# Line 700  class Symbol(object): Line 779  class Symbol(object):
779         """         """
780         returns -self.         returns -self.
781    
782         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
783         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
784         """         """
785         return self*(-1.)         return self*(-1.)
# Line 709  class Symbol(object): Line 788  class Symbol(object):
788         """         """
789         returns +self.         returns +self.
790    
791         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
792         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
793         """         """
794         return self*(1.)         return self*(1.)
795    
796     def __abs__(self):     def __abs__(self):
797         """         """
798         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
799         """         """
800         return Abs_Symbol(self)         return Abs_Symbol(self)
801    
# Line 726  class Symbol(object): Line 805  class Symbol(object):
805    
806         @param other: object to be added to this object         @param other: object to be added to this object
807         @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}.
808         @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}
809         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
810         """         """
811         return add(self,other)         return add(self,other)
# Line 737  class Symbol(object): Line 816  class Symbol(object):
816    
817         @param other: object this object is added to         @param other: object this object is added to
818         @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}.
819         @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
820         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
821         """         """
822         return add(other,self)         return add(other,self)
# Line 748  class Symbol(object): Line 827  class Symbol(object):
827    
828         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
829         @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}.
830         @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
831         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
832         """         """
833         return add(self,-other)         return add(self,-other)
# Line 759  class Symbol(object): Line 838  class Symbol(object):
838    
839         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
840         @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}.
841         @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}.
842         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
843         """         """
844         return add(-self,other)         return add(-self,other)
# Line 770  class Symbol(object): Line 849  class Symbol(object):
849    
850         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
851         @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}.
852         @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}.
853         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
854         """         """
855         return mult(self,other)         return mult(self,other)
# Line 781  class Symbol(object): Line 860  class Symbol(object):
860    
861         @param other: object this object is multiplied with         @param other: object this object is multiplied with
862         @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}.
863         @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.
864         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
865         """         """
866         return mult(other,self)         return mult(other,self)
# Line 792  class Symbol(object): Line 871  class Symbol(object):
871    
872         @param other: object dividing this object         @param other: object dividing this object
873         @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}.
874         @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}
875         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
876         """         """
877         return quotient(self,other)         return quotient(self,other)
# Line 803  class Symbol(object): Line 882  class Symbol(object):
882    
883         @param other: object dividing this object         @param other: object dividing this object
884         @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}.
885         @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
886         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
887         """         """
888         return quotient(other,self)         return quotient(other,self)
# Line 814  class Symbol(object): Line 893  class Symbol(object):
893    
894         @param other: exponent         @param other: exponent
895         @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}.
896         @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}
897         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
898         """         """
899         return power(self,other)         return power(self,other)
# Line 825  class Symbol(object): Line 904  class Symbol(object):
904    
905         @param other: basis         @param other: basis
906         @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}.
907         @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
908         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
909         """         """
910         return power(other,self)         return power(other,self)
911    
912       def __getitem__(self,index):
913           """
914           returns the slice defined by index
915    
916           @param index: defines a
917           @type index: C{slice} or C{int} or a C{tuple} of them
918           @return: a L{Symbol} representing the slice defined by index
919           @rtype: L{DependendSymbol}
920           """
921           return GetSlice_Symbol(self,index)
922    
923  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
924     """     """
925     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.
926     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  
927        
928     Example:     Example::
929        
930     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
931     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
932     print u1==u2       print u1==u2
933     False       False
934        
935        but       but::
936    
937     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
938     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
939     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
940     print u1==u2, u1==u3       print u1==u2, u1==u3
941     True False       True False
942    
943     @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.
944     """     """
# Line 880  class DependendSymbol(Symbol): Line 970  class DependendSymbol(Symbol):
970  #=========================================================  #=========================================================
971  #  Unary operations prserving the shape  #  Unary operations prserving the shape
972  #========================================================  #========================================================
973    class GetSlice_Symbol(DependendSymbol):
974       """
975       L{Symbol} representing getting a slice for a L{Symbol}
976       """
977       def __init__(self,arg,index):
978          """
979          initialization of wherePositive L{Symbol} with argument arg
980          @param arg: argument
981          @type arg: L{Symbol}.
982          @param index: defines index
983          @type index: C{slice} or C{int} or a C{tuple} of them
984          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
985          @raises ValueError: if a step is given
986          """
987          if not isinstance(index,tuple): index=(index,)
988          if len(index)>arg.getRank():
989               raise IndexError,"GetSlice_Symbol: index out of range."
990          sh=()
991          index2=()
992          for i in range(len(index)):
993             ix=index[i]
994             if isinstance(ix,int):
995                if ix<0 or ix>=arg.getShape()[i]:
996                   raise ValueError,"GetSlice_Symbol: index out of range."
997                index2=index2+(ix,)
998             else:
999               if not ix.step==None:
1000                 raise ValueError,"GetSlice_Symbol: steping is not supported."
1001               if ix.start==None:
1002                  s=0
1003               else:
1004                  s=ix.start
1005               if ix.stop==None:
1006                  e=arg.getShape()[i]
1007               else:
1008                  e=ix.stop
1009                  if e>arg.getShape()[i]:
1010                     raise IndexError,"GetSlice_Symbol: index out of range."
1011               index2=index2+(slice(s,e),)
1012               if e>s:
1013                   sh=sh+(e-s,)
1014               elif s>e:
1015                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
1016          for i in range(len(index),arg.getRank()):
1017              index2=index2+(slice(0,arg.getShape()[i]),)
1018              sh=sh+(arg.getShape()[i],)
1019          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
1020    
1021       def getMyCode(self,argstrs,format="escript"):
1022          """
1023          returns a program code that can be used to evaluate the symbol.
1024    
1025          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1026          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1027          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1028          @type format: C{str}
1029          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1030          @rtype: C{str}
1031          @raise NotImplementedError: if the requested format is not available
1032          """
1033          if format=="escript" or format=="str"  or format=="text":
1034             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
1035          else:
1036             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
1037    
1038       def substitute(self,argvals):
1039          """
1040          assigns new values to symbols in the definition of the symbol.
1041          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1042    
1043          @param argvals: new values assigned to symbols
1044          @type argvals: C{dict} with keywords of type L{Symbol}.
1045          @return: result of the substitution process. Operations are executed as much as possible.
1046          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1047          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1048          """
1049          if argvals.has_key(self):
1050             arg=argvals[self]
1051             if self.isAppropriateValue(arg):
1052                return arg
1053             else:
1054                raise TypeError,"%s: new value is not appropriate."%str(self)
1055          else:
1056             args=self.getSubstitutedArguments(argvals)
1057             arg=args[0]
1058             index=args[1]
1059             return arg.__getitem__(index)
1060    
1061  def log10(arg):  def log10(arg):
1062     """     """
1063     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
1064    
1065     @param arg: argument     @param arg: argument
1066     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1067     @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.
1068     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1069     """     """
1070     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 908  def wherePositive(arg): Line 1086  def wherePositive(arg):
1086    
1087     @param arg: argument     @param arg: argument
1088     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1089     @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.
1090     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1091     """     """
1092     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1093        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1094        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1095        return out        return out
1096     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1097        return arg._wherePositive()        return arg._wherePositive()
# Line 954  class WherePositive_Symbol(DependendSymb Line 1132  class WherePositive_Symbol(DependendSymb
1132        @type format: C{str}        @type format: C{str}
1133        @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.
1134        @rtype: C{str}        @rtype: C{str}
1135        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1136        """        """
1137        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1138            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 990  def whereNegative(arg): Line 1168  def whereNegative(arg):
1168    
1169     @param arg: argument     @param arg: argument
1170     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1171     @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.
1172     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1173     """     """
1174     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1175        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1176        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1177        return out        return out
1178     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1179        return arg._whereNegative()        return arg._whereNegative()
# Line 1036  class WhereNegative_Symbol(DependendSymb Line 1214  class WhereNegative_Symbol(DependendSymb
1214        @type format: C{str}        @type format: C{str}
1215        @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.
1216        @rtype: C{str}        @rtype: C{str}
1217        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1218        """        """
1219        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1220            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1072  def whereNonNegative(arg): Line 1250  def whereNonNegative(arg):
1250    
1251     @param arg: argument     @param arg: argument
1252     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1253     @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.
1254     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1255     """     """
1256     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1257        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1258        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1259        return out        return out
1260     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1261        return arg._whereNonNegative()        return arg._whereNonNegative()
# Line 1102  def whereNonPositive(arg): Line 1280  def whereNonPositive(arg):
1280    
1281     @param arg: argument     @param arg: argument
1282     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1283     @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.
1284     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1285     """     """
1286     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1287        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1288        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1289        return out        return out
1290     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1291        return arg._whereNonPositive()        return arg._whereNonPositive()
# Line 1134  def whereZero(arg,tol=0.): Line 1312  def whereZero(arg,tol=0.):
1312     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1313     @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.
1314     @type tol: C{float}     @type tol: C{float}
1315     @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.
1316     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1317     """     """
1318     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1319        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1320        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1321        return out        return out
1322     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1323        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1324     elif isinstance(arg,float):     elif isinstance(arg,float):
1325        if abs(arg)<=tol:        if abs(arg)<=tol:
1326          return 1.          return 1.
# Line 1183  class WhereZero_Symbol(DependendSymbol): Line 1358  class WhereZero_Symbol(DependendSymbol):
1358        @type format: C{str}        @type format: C{str}
1359        @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.
1360        @rtype: C{str}        @rtype: C{str}
1361        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1362        """        """
1363        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1364           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1217  def whereNonZero(arg,tol=0.): Line 1392  def whereNonZero(arg,tol=0.):
1392    
1393     @param arg: argument     @param arg: argument
1394     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1395     @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.
1396     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1397     """     """
1398     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1399        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1400        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1401        return out        return out
1402     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1403        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1404     elif isinstance(arg,float):     elif isinstance(arg,float):
1405        if abs(arg)>tol:        if abs(arg)>tol:
1406          return 1.          return 1.
# Line 1244  def whereNonZero(arg,tol=0.): Line 1416  def whereNonZero(arg,tol=0.):
1416     else:     else:
1417        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1418    
1419    def erf(arg):
1420       """
1421       returns erf of argument arg
1422    
1423       @param arg: argument
1424       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1425       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1426       @raises TypeError: if the type of the argument is not expected.
1427       """
1428       if isinstance(arg,escript.Data):
1429          return arg._erf()
1430       else:
1431          raise TypeError,"erf: Unknown argument type."
1432    
1433  def sin(arg):  def sin(arg):
1434     """     """
1435     returns sine of argument arg     returns sine of argument arg
1436    
1437     @param arg: argument     @param arg: argument
1438     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1439     @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.
1440     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1441     """     """
1442     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1288  class Sin_Symbol(DependendSymbol): Line 1474  class Sin_Symbol(DependendSymbol):
1474        @type format: C{str}        @type format: C{str}
1475        @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.
1476        @rtype: C{str}        @rtype: C{str}
1477        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1478        """        """
1479        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1480            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1340  def cos(arg): Line 1526  def cos(arg):
1526    
1527     @param arg: argument     @param arg: argument
1528     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1529     @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.
1530     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1531     """     """
1532     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1378  class Cos_Symbol(DependendSymbol): Line 1564  class Cos_Symbol(DependendSymbol):
1564        @type format: C{str}        @type format: C{str}
1565        @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.
1566        @rtype: C{str}        @rtype: C{str}
1567        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1568        """        """
1569        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1570            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1430  def tan(arg): Line 1616  def tan(arg):
1616    
1617     @param arg: argument     @param arg: argument
1618     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1619     @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.
1620     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1621     """     """
1622     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1468  class Tan_Symbol(DependendSymbol): Line 1654  class Tan_Symbol(DependendSymbol):
1654        @type format: C{str}        @type format: C{str}
1655        @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.
1656        @rtype: C{str}        @rtype: C{str}
1657        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1658        """        """
1659        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1660            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1520  def asin(arg): Line 1706  def asin(arg):
1706    
1707     @param arg: argument     @param arg: argument
1708     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1709     @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.
1710     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1711     """     """
1712     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1558  class Asin_Symbol(DependendSymbol): Line 1744  class Asin_Symbol(DependendSymbol):
1744        @type format: C{str}        @type format: C{str}
1745        @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.
1746        @rtype: C{str}        @rtype: C{str}
1747        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1748        """        """
1749        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1750            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1610  def acos(arg): Line 1796  def acos(arg):
1796    
1797     @param arg: argument     @param arg: argument
1798     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1799     @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.
1800     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1801     """     """
1802     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1648  class Acos_Symbol(DependendSymbol): Line 1834  class Acos_Symbol(DependendSymbol):
1834        @type format: C{str}        @type format: C{str}
1835        @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.
1836        @rtype: C{str}        @rtype: C{str}
1837        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1838        """        """
1839        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1840            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1700  def atan(arg): Line 1886  def atan(arg):
1886    
1887     @param arg: argument     @param arg: argument
1888     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1889     @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.
1890     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1891     """     """
1892     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1738  class Atan_Symbol(DependendSymbol): Line 1924  class Atan_Symbol(DependendSymbol):
1924        @type format: C{str}        @type format: C{str}
1925        @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.
1926        @rtype: C{str}        @rtype: C{str}
1927        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1928        """        """
1929        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1930            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1790  def sinh(arg): Line 1976  def sinh(arg):
1976    
1977     @param arg: argument     @param arg: argument
1978     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1979     @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.
1980     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1981     """     """
1982     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1828  class Sinh_Symbol(DependendSymbol): Line 2014  class Sinh_Symbol(DependendSymbol):
2014        @type format: C{str}        @type format: C{str}
2015        @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.
2016        @rtype: C{str}        @rtype: C{str}
2017        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2018        """        """
2019        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2020            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1880  def cosh(arg): Line 2066  def cosh(arg):
2066    
2067     @param arg: argument     @param arg: argument
2068     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2069     @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.
2070     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2071     """     """
2072     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1918  class Cosh_Symbol(DependendSymbol): Line 2104  class Cosh_Symbol(DependendSymbol):
2104        @type format: C{str}        @type format: C{str}
2105        @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.
2106        @rtype: C{str}        @rtype: C{str}
2107        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2108        """        """
2109        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2110            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1970  def tanh(arg): Line 2156  def tanh(arg):
2156    
2157     @param arg: argument     @param arg: argument
2158     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2159     @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.
2160     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2161     """     """
2162     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2008  class Tanh_Symbol(DependendSymbol): Line 2194  class Tanh_Symbol(DependendSymbol):
2194        @type format: C{str}        @type format: C{str}
2195        @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.
2196        @rtype: C{str}        @rtype: C{str}
2197        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2198        """        """
2199        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2200            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2060  def asinh(arg): Line 2246  def asinh(arg):
2246    
2247     @param arg: argument     @param arg: argument
2248     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2249     @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.
2250     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2251     """     """
2252     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2098  class Asinh_Symbol(DependendSymbol): Line 2284  class Asinh_Symbol(DependendSymbol):
2284        @type format: C{str}        @type format: C{str}
2285        @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.
2286        @rtype: C{str}        @rtype: C{str}
2287        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2288        """        """
2289        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2290            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2150  def acosh(arg): Line 2336  def acosh(arg):
2336    
2337     @param arg: argument     @param arg: argument
2338     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2339     @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.
2340     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2341     """     """
2342     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2188  class Acosh_Symbol(DependendSymbol): Line 2374  class Acosh_Symbol(DependendSymbol):
2374        @type format: C{str}        @type format: C{str}
2375        @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.
2376        @rtype: C{str}        @rtype: C{str}
2377        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2378        """        """
2379        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2380            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2240  def atanh(arg): Line 2426  def atanh(arg):
2426    
2427     @param arg: argument     @param arg: argument
2428     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2429     @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.
2430     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2431     """     """
2432     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2278  class Atanh_Symbol(DependendSymbol): Line 2464  class Atanh_Symbol(DependendSymbol):
2464        @type format: C{str}        @type format: C{str}
2465        @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.
2466        @rtype: C{str}        @rtype: C{str}
2467        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2468        """        """
2469        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2470            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2330  def exp(arg): Line 2516  def exp(arg):
2516    
2517     @param arg: argument     @param arg: argument
2518     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2519     @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.
2520     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2521     """     """
2522     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2368  class Exp_Symbol(DependendSymbol): Line 2554  class Exp_Symbol(DependendSymbol):
2554        @type format: C{str}        @type format: C{str}
2555        @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.
2556        @rtype: C{str}        @rtype: C{str}
2557        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2558        """        """
2559        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2560            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2420  def sqrt(arg): Line 2606  def sqrt(arg):
2606    
2607     @param arg: argument     @param arg: argument
2608     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2609     @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.
2610     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2611     """     """
2612     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2458  class Sqrt_Symbol(DependendSymbol): Line 2644  class Sqrt_Symbol(DependendSymbol):
2644        @type format: C{str}        @type format: C{str}
2645        @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.
2646        @rtype: C{str}        @rtype: C{str}
2647        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2648        """        """
2649        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2650            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2510  def log(arg): Line 2696  def log(arg):
2696    
2697     @param arg: argument     @param arg: argument
2698     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2699     @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.
2700     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2701     """     """
2702     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2548  class Log_Symbol(DependendSymbol): Line 2734  class Log_Symbol(DependendSymbol):
2734        @type format: C{str}        @type format: C{str}
2735        @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.
2736        @rtype: C{str}        @rtype: C{str}
2737        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2738        """        """
2739        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2740            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2600  def sign(arg): Line 2786  def sign(arg):
2786    
2787     @param arg: argument     @param arg: argument
2788     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2789     @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.
2790     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2791     """     """
2792     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2648  class Abs_Symbol(DependendSymbol): Line 2834  class Abs_Symbol(DependendSymbol):
2834        @type format: C{str}        @type format: C{str}
2835        @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.
2836        @rtype: C{str}        @rtype: C{str}
2837        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2838        """        """
2839        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2840            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2700  def minval(arg): Line 2886  def minval(arg):
2886    
2887     @param arg: argument     @param arg: argument
2888     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2889     @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.
2890     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2891     """     """
2892     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2741  class Minval_Symbol(DependendSymbol): Line 2927  class Minval_Symbol(DependendSymbol):
2927        @type format: C{str}        @type format: C{str}
2928        @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.
2929        @rtype: C{str}        @rtype: C{str}
2930        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2931        """        """
2932        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2933            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2777  def maxval(arg): Line 2963  def maxval(arg):
2963    
2964     @param arg: argument     @param arg: argument
2965     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2966     @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.
2967     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2968     """     """
2969     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2818  class Maxval_Symbol(DependendSymbol): Line 3004  class Maxval_Symbol(DependendSymbol):
3004        @type format: C{str}        @type format: C{str}
3005        @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.
3006        @rtype: C{str}        @rtype: C{str}
3007        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3008        """        """
3009        if isinstance(argstrs,list):        if isinstance(argstrs,list):
3010            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2854  def length(arg): Line 3040  def length(arg):
3040    
3041     @param arg: argument     @param arg: argument
3042     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3043     @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.
3044     """     """
3045     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
3046    
# Line 2864  def trace(arg,axis_offset=0): Line 3050  def trace(arg,axis_offset=0):
3050    
3051     @param arg: argument     @param arg: argument
3052     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3053     @param axis_offset: 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     @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
3054                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3055     @type axis_offset: C{int}     @type axis_offset: C{int}
3056     @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.     @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
3057     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
# Line 2873  def trace(arg,axis_offset=0): Line 3059  def trace(arg,axis_offset=0):
3059     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
3060        sh=arg.shape        sh=arg.shape
3061        if len(sh)<2:        if len(sh)<2:
3062          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3063        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
3064          raise ValueError,"trace: axis_offset must be between 0 and %s"%len(sh)-2          raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
3065        s1=1        s1=1
3066        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
3067        s2=1        s2=1
3068        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3069        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
3070          raise ValueError,"trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3071        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
3072        out=numarray.zeros([s1,s2],numarray.Float)        out=numarray.zeros([s1,s2],numarray.Float64)
3073        for i1 in range(s1):        for i1 in range(s1):
3074          for i2 in range(s2):          for i2 in range(s2):
3075              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
3076        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3077        return out        return out
3078     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3079        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
3080            raise ValueError,"rank of argument must be greater than 1"
3081          if axis_offset<0 or axis_offset>arg.getRank()-2:
3082            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3083          s=list(arg.getShape())        
3084          if not s[axis_offset] == s[axis_offset+1]:
3085            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3086          return arg._trace(axis_offset)
3087     elif isinstance(arg,float):     elif isinstance(arg,float):
3088        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
3089     elif isinstance(arg,int):     elif isinstance(arg,int):
3090        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
3091     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3092        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
3093     else:     else:
3094        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
3095    
 def escript_trace(arg,axis_offset): # this should be escript._trace  
       "arg si a Data objects!!!"  
       if arg.getRank()<2:  
         raise ValueError,"escript_trace: rank of argument must be greater than 1"  
       if axis_offset<0 or axis_offset>arg.getRank()-2:  
         raise ValueError,"escript_trace: axis_offset must be between 0 and %s"%arg.getRank()-2  
       s=list(arg.getShape())          
       if not s[axis_offset] == s[axis_offset+1]:  
         raise ValueError,"escript_trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)  
       out=escript.Data(0.,tuple(s[0:axis_offset]+s[axis_offset+2:]),arg.getFunctionSpace())  
       if arg.getRank()==2:  
          for i0 in range(s[0]):  
             out+=arg[i0,i0]  
       elif arg.getRank()==3:  
          if axis_offset==0:  
             for i0 in range(s[0]):  
                   for i2 in range(s[2]):  
                          out[i2]+=arg[i0,i0,i2]  
          elif axis_offset==1:  
             for i0 in range(s[0]):  
                for i1 in range(s[1]):  
                          out[i0]+=arg[i0,i1,i1]  
       elif arg.getRank()==4:  
          if axis_offset==0:  
             for i0 in range(s[0]):  
                   for i2 in range(s[2]):  
                      for i3 in range(s[3]):  
                          out[i2,i3]+=arg[i0,i0,i2,i3]  
          elif axis_offset==1:  
             for i0 in range(s[0]):  
                for i1 in range(s[1]):  
                      for i3 in range(s[3]):  
                          out[i0,i3]+=arg[i0,i1,i1,i3]  
          elif axis_offset==2:  
             for i0 in range(s[0]):  
                for i1 in range(s[1]):  
                   for i2 in range(s[2]):  
                          out[i0,i1]+=arg[i0,i1,i2,i2]  
       return out  
3096  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
3097     """     """
3098     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 2948  class Trace_Symbol(DependendSymbol): Line 3102  class Trace_Symbol(DependendSymbol):
3102        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3103        @param arg: argument of function        @param arg: argument of function
3104        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3105        @param axis_offset: 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        @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
3106                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3107        @type axis_offset: C{int}        @type axis_offset: C{int}
3108        """        """
3109        if arg.getRank()<2:        if arg.getRank()<2:
3110          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3111        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3112          raise ValueError,"Trace_Symbol: axis_offset must be between 0 and %s"%arg.getRank()-2          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3113        s=list(arg.getShape())                s=list(arg.getShape())        
3114        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3115          raise ValueError,"Trace_Symbol: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3116        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3117    
3118     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 2971  class Trace_Symbol(DependendSymbol): Line 3125  class Trace_Symbol(DependendSymbol):
3125        @type format: C{str}        @type format: C{str}
3126        @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.
3127        @rtype: C{str}        @rtype: C{str}
3128        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3129        """        """
3130        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3131           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3013  class Trace_Symbol(DependendSymbol): Line 3167  class Trace_Symbol(DependendSymbol):
3167        else:        else:
3168           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3169    
3170    def transpose(arg,axis_offset=None):
3171       """
3172       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3173    
3174       @param arg: argument
3175       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3176       @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.
3177                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3178       @type axis_offset: C{int}
3179       @return: transpose of arg
3180       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3181       """
3182       if isinstance(arg,numarray.NumArray):
3183          if axis_offset==None: axis_offset=int(arg.rank/2)
3184          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3185       elif isinstance(arg,escript.Data):
3186          r=arg.getRank()
3187          if axis_offset==None: axis_offset=int(r/2)
3188          if axis_offset<0 or axis_offset>r:
3189            raise ValueError,"axis_offset must be between 0 and %s"%r
3190          return arg._transpose(axis_offset)
3191       elif isinstance(arg,float):
3192          if not ( axis_offset==0 or axis_offset==None):
3193            raise ValueError,"axis_offset must be 0 for float argument"
3194          return arg
3195       elif isinstance(arg,int):
3196          if not ( axis_offset==0 or axis_offset==None):
3197            raise ValueError,"axis_offset must be 0 for int argument"
3198          return float(arg)
3199       elif isinstance(arg,Symbol):
3200          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3201          return Transpose_Symbol(arg,axis_offset)
3202       else:
3203          raise TypeError,"Unknown argument type."
3204    
3205    class Transpose_Symbol(DependendSymbol):
3206       """
3207       L{Symbol} representing the result of the transpose function
3208       """
3209       def __init__(self,arg,axis_offset=None):
3210          """
3211          initialization of transpose L{Symbol} with argument arg
3212    
3213          @param arg: argument of function
3214          @type arg: L{Symbol}.
3215          @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.
3216                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3217          @type axis_offset: C{int}
3218          """
3219          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3220          if axis_offset<0 or axis_offset>arg.getRank():
3221            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3222          s=arg.getShape()
3223          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3224    
3225       def getMyCode(self,argstrs,format="escript"):
3226          """
3227          returns a program code that can be used to evaluate the symbol.
3228    
3229          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3230          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3231          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3232          @type format: C{str}
3233          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3234          @rtype: C{str}
3235          @raise NotImplementedError: if the requested format is not available
3236          """
3237          if format=="escript" or format=="str"  or format=="text":
3238             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3239          else:
3240             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3241    
3242       def substitute(self,argvals):
3243          """
3244          assigns new values to symbols in the definition of the symbol.
3245          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3246    
3247          @param argvals: new values assigned to symbols
3248          @type argvals: C{dict} with keywords of type L{Symbol}.
3249          @return: result of the substitution process. Operations are executed as much as possible.
3250          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3251          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3252          """
3253          if argvals.has_key(self):
3254             arg=argvals[self]
3255             if self.isAppropriateValue(arg):
3256                return arg
3257             else:
3258                raise TypeError,"%s: new value is not appropriate."%str(self)
3259          else:
3260             arg=self.getSubstitutedArguments(argvals)
3261             return transpose(arg[0],axis_offset=arg[1])
3262    
3263       def diff(self,arg):
3264          """
3265          differential of this object
3266    
3267          @param arg: the derivative is calculated with respect to arg
3268          @type arg: L{escript.Symbol}
3269          @return: derivative with respect to C{arg}
3270          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3271          """
3272          if arg==self:
3273             return identity(self.getShape())
3274          else:
3275             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3276    
3277    def swap_axes(arg,axis0=0,axis1=1):
3278       """
3279       returns the swap of arg by swaping the components axis0 and axis1
3280    
3281       @param arg: argument
3282       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3283       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3284       @type axis0: C{int}
3285       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3286       @type axis1: C{int}
3287       @return: C{arg} with swaped components
3288       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3289       """
3290       if axis0 > axis1:
3291          axis0,axis1=axis1,axis0
3292       if isinstance(arg,numarray.NumArray):
3293          return numarray.swapaxes(arg,axis0,axis1)
3294       elif isinstance(arg,escript.Data):
3295          return arg._swap_axes(axis0,axis1)
3296       elif isinstance(arg,float):
3297          raise TyepError,"float argument is not supported."
3298       elif isinstance(arg,int):
3299          raise TyepError,"int argument is not supported."
3300       elif isinstance(arg,Symbol):
3301          return SwapAxes_Symbol(arg,axis0,axis1)
3302       else:
3303          raise TypeError,"Unknown argument type."
3304    
3305    class SwapAxes_Symbol(DependendSymbol):
3306       """
3307       L{Symbol} representing the result of the swap function
3308       """
3309       def __init__(self,arg,axis0=0,axis1=1):
3310          """
3311          initialization of swap L{Symbol} with argument arg
3312    
3313          @param arg: argument
3314          @type arg: L{Symbol}.
3315          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3316          @type axis0: C{int}
3317          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3318          @type axis1: C{int}
3319          """
3320          if arg.getRank()<2:
3321             raise ValueError,"argument must have at least rank 2."
3322          if axis0<0 or axis0>arg.getRank()-1:
3323             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3324          if axis1<0 or axis1>arg.getRank()-1:
3325             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3326          if axis0 == axis1:
3327             raise ValueError,"axis indices must be different."
3328          if axis0 > axis1:
3329             axis0,axis1=axis1,axis0
3330          s=arg.getShape()
3331          s_out=[]
3332          for i in range(len(s)):
3333             if i == axis0:
3334                s_out.append(s[axis1])
3335             elif i == axis1:
3336                s_out.append(s[axis0])
3337             else:
3338                s_out.append(s[i])
3339          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3340    
3341       def getMyCode(self,argstrs,format="escript"):
3342          """
3343          returns a program code that can be used to evaluate the symbol.
3344    
3345          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3346          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3347          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3348          @type format: C{str}
3349          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3350          @rtype: C{str}
3351          @raise NotImplementedError: if the requested format is not available
3352          """
3353          if format=="escript" or format=="str"  or format=="text":
3354             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3355          else:
3356             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3357    
3358       def substitute(self,argvals):
3359          """
3360          assigns new values to symbols in the definition of the symbol.
3361          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3362    
3363          @param argvals: new values assigned to symbols
3364          @type argvals: C{dict} with keywords of type L{Symbol}.
3365          @return: result of the substitution process. Operations are executed as much as possible.
3366          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3367          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3368          """
3369          if argvals.has_key(self):
3370             arg=argvals[self]
3371             if self.isAppropriateValue(arg):
3372                return arg
3373             else:
3374                raise TypeError,"%s: new value is not appropriate."%str(self)
3375          else:
3376             arg=self.getSubstitutedArguments(argvals)
3377             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3378    
3379       def diff(self,arg):
3380          """
3381          differential of this object
3382    
3383          @param arg: the derivative is calculated with respect to arg
3384          @type arg: L{escript.Symbol}
3385          @return: derivative with respect to C{arg}
3386          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3387          """
3388          if arg==self:
3389             return identity(self.getShape())
3390          else:
3391             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3392    
3393    def symmetric(arg):
3394        """
3395        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3396    
3397        @param arg: square matrix. Must have rank 2 or 4 and be square.
3398        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3399        @return: symmetric part of arg
3400        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3401        """
3402        if isinstance(arg,numarray.NumArray):
3403          if arg.rank==2:
3404            if not (arg.shape[0]==arg.shape[1]):
3405               raise ValueError,"argument must be square."
3406          elif arg.rank==4:
3407            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3408               raise ValueError,"argument must be square."
3409          else:
3410            raise ValueError,"rank 2 or 4 is required."
3411          return (arg+transpose(arg))/2
3412        elif isinstance(arg,escript.Data):
3413          if arg.getRank()==2:
3414            if not (arg.getShape()[0]==arg.getShape()[1]):
3415               raise ValueError,"argument must be square."
3416            return arg._symmetric()
3417          elif arg.getRank()==4:
3418            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3419               raise ValueError,"argument must be square."
3420            return arg._symmetric()
3421          else:
3422            raise ValueError,"rank 2 or 4 is required."
3423        elif isinstance(arg,float):
3424          return arg
3425        elif isinstance(arg,int):
3426          return float(arg)
3427        elif isinstance(arg,Symbol):
3428          if arg.getRank()==2:
3429            if not (arg.getShape()[0]==arg.getShape()[1]):
3430               raise ValueError,"argument must be square."
3431          elif arg.getRank()==4:
3432            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3433               raise ValueError,"argument must be square."
3434          else:
3435            raise ValueError,"rank 2 or 4 is required."
3436          return (arg+transpose(arg))/2
3437        else:
3438          raise TypeError,"symmetric: Unknown argument type."
3439    
3440    def nonsymmetric(arg):
3441        """
3442        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3443    
3444        @param arg: square matrix. Must have rank 2 or 4 and be square.
3445        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3446        @return: nonsymmetric part of arg
3447        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3448        """
3449        if isinstance(arg,numarray.NumArray):
3450          if arg.rank==2:
3451            if not (arg.shape[0]==arg.shape[1]):
3452               raise ValueError,"nonsymmetric: argument must be square."
3453          elif arg.rank==4:
3454            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3455               raise ValueError,"nonsymmetric: argument must be square."
3456          else:
3457            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3458          return (arg-transpose(arg))/2
3459        elif isinstance(arg,escript.Data):
3460          if arg.getRank()==2:
3461            if not (arg.getShape()[0]==arg.getShape()[1]):
3462               raise ValueError,"argument must be square."
3463            return arg._nonsymmetric()
3464          elif arg.getRank()==4:
3465            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3466               raise ValueError,"argument must be square."
3467            return arg._nonsymmetric()
3468          else:
3469            raise ValueError,"rank 2 or 4 is required."
3470        elif isinstance(arg,float):
3471          return arg
3472        elif isinstance(arg,int):
3473          return float(arg)
3474        elif isinstance(arg,Symbol):
3475          if arg.getRank()==2:
3476            if not (arg.getShape()[0]==arg.getShape()[1]):
3477               raise ValueError,"nonsymmetric: argument must be square."
3478          elif arg.getRank()==4:
3479            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3480               raise ValueError,"nonsymmetric: argument must be square."
3481          else:
3482            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3483          return (arg-transpose(arg))/2
3484        else:
3485          raise TypeError,"nonsymmetric: Unknown argument type."
3486    
3487    def inverse(arg):
3488        """
3489        returns the inverse of the square matrix arg.
3490    
3491        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3492        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3493        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3494        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3495        @note: for L{escript.Data} objects the dimension is restricted to 3.
3496        """
3497        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3498        if isinstance(arg,numarray.NumArray):
3499          return numarray.linear_algebra.inverse(arg)
3500        elif isinstance(arg,escript.Data):
3501          return escript_inverse(arg)
3502        elif isinstance(arg,float):
3503          return 1./arg
3504        elif isinstance(arg,int):
3505          return 1./float(arg)
3506        elif isinstance(arg,Symbol):
3507          return Inverse_Symbol(arg)
3508        else:
3509          raise TypeError,"inverse: Unknown argument type."
3510    
3511    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3512          "arg is a Data objects!!!"
3513          if not arg.getRank()==2:
3514            raise ValueError,"escript_inverse: argument must have rank 2"
3515          s=arg.getShape()      
3516          if not s[0] == s[1]:
3517            raise ValueError,"escript_inverse: argument must be a square matrix."
3518          out=escript.Data(0.,s,arg.getFunctionSpace())
3519          if s[0]==1:
3520              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3521                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3522              out[0,0]=1./arg[0,0]
3523          elif s[0]==2:
3524              A11=arg[0,0]
3525              A12=arg[0,1]
3526              A21=arg[1,0]
3527              A22=arg[1,1]
3528              D = A11*A22-A12*A21
3529              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3530                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3531              D=1./D
3532              out[0,0]= A22*D
3533              out[1,0]=-A21*D
3534              out[0,1]=-A12*D
3535              out[1,1]= A11*D
3536          elif s[0]==3:
3537              A11=arg[0,0]
3538              A21=arg[1,0]
3539              A31=arg[2,0]
3540              A12=arg[0,1]
3541              A22=arg[1,1]
3542              A32=arg[2,1]
3543              A13=arg[0,2]
3544              A23=arg[1,2]
3545              A33=arg[2,2]
3546              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3547              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3548                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3549              D=1./D
3550              out[0,0]=(A22*A33-A23*A32)*D
3551              out[1,0]=(A31*A23-A21*A33)*D
3552              out[2,0]=(A21*A32-A31*A22)*D
3553              out[0,1]=(A13*A32-A12*A33)*D
3554              out[1,1]=(A11*A33-A31*A13)*D
3555              out[2,1]=(A12*A31-A11*A32)*D
3556              out[0,2]=(A12*A23-A13*A22)*D
3557              out[1,2]=(A13*A21-A11*A23)*D
3558              out[2,2]=(A11*A22-A12*A21)*D
3559          else:
3560             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3561          return out
3562    
3563    class Inverse_Symbol(DependendSymbol):
3564       """
3565       L{Symbol} representing the result of the inverse function
3566       """
3567       def __init__(self,arg):
3568          """
3569          initialization of inverse L{Symbol} with argument arg
3570          @param arg: argument of function
3571          @type arg: L{Symbol}.
3572          """
3573          if not arg.getRank()==2:
3574            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3575          s=arg.getShape()
3576          if not s[0] == s[1]:
3577            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3578          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3579    
3580       def getMyCode(self,argstrs,format="escript"):
3581          """
3582          returns a program code that can be used to evaluate the symbol.
3583    
3584          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3585          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3586          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3587          @type format: C{str}
3588          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3589          @rtype: C{str}
3590          @raise NotImplementedError: if the requested format is not available
3591          """
3592          if format=="escript" or format=="str"  or format=="text":
3593             return "inverse(%s)"%argstrs[0]
3594          else:
3595             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3596    
3597       def substitute(self,argvals):
3598          """
3599          assigns new values to symbols in the definition of the symbol.
3600          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3601    
3602          @param argvals: new values assigned to symbols
3603          @type argvals: C{dict} with keywords of type L{Symbol}.
3604          @return: result of the substitution process. Operations are executed as much as possible.
3605          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3606          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3607          """
3608          if argvals.has_key(self):
3609             arg=argvals[self]
3610             if self.isAppropriateValue(arg):
3611                return arg
3612             else:
3613                raise TypeError,"%s: new value is not appropriate."%str(self)
3614          else:
3615             arg=self.getSubstitutedArguments(argvals)
3616             return inverse(arg[0])
3617    
3618       def diff(self,arg):
3619          """
3620          differential of this object
3621    
3622          @param arg: the derivative is calculated with respect to arg
3623          @type arg: L{escript.Symbol}
3624          @return: derivative with respect to C{arg}
3625          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3626          """
3627          if arg==self:
3628             return identity(self.getShape())
3629          else:
3630             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3631    
3632    def eigenvalues(arg):
3633        """
3634        returns the eigenvalues of the square matrix arg.
3635    
3636        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3637                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3638        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3639        @return: the eigenvalues in increasing order.
3640        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3641        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3642        """
3643        if isinstance(arg,numarray.NumArray):
3644          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3645          out.sort()
3646          return out
3647        elif isinstance(arg,escript.Data):
3648          return arg._eigenvalues()
3649        elif isinstance(arg,Symbol):
3650          if not arg.getRank()==2:
3651            raise ValueError,"eigenvalues: argument must have rank 2"
3652          s=arg.getShape()      
3653          if not s[0] == s[1]:
3654            raise ValueError,"eigenvalues: argument must be a square matrix."
3655          if s[0]==1:
3656              return arg[0]
3657          elif s[0]==2:
3658              arg1=symmetric(arg)
3659              A11=arg1[0,0]
3660              A12=arg1[0,1]
3661              A22=arg1[1,1]
3662              trA=(A11+A22)/2.
3663              A11-=trA
3664              A22-=trA
3665              s=sqrt(A12**2-A11*A22)
3666              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3667          elif s[0]==3:
3668              arg1=symmetric(arg)
3669              A11=arg1[0,0]
3670              A12=arg1[0,1]
3671              A22=arg1[1,1]
3672              A13=arg1[0,2]
3673              A23=arg1[1,2]
3674              A33=arg1[2,2]
3675              trA=(A11+A22+A33)/3.
3676              A11-=trA
3677              A22-=trA
3678              A33-=trA
3679              A13_2=A13**2
3680              A23_2=A23**2
3681              A12_2=A12**2
3682              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3683              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3684              sq_p=sqrt(p/3.)
3685              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
3686              sq_p*=2.
3687              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3688               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3689               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3690              return trA+sq_p*f
3691          else:
3692             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3693        elif isinstance(arg,float):
3694          return arg
3695        elif isinstance(arg,int):
3696          return float(arg)
3697        else:
3698          raise TypeError,"eigenvalues: Unknown argument type."
3699    
3700    def eigenvalues_and_eigenvectors(arg):
3701        """
3702        returns the eigenvalues and eigenvectors of the square matrix arg.
3703    
3704        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3705                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3706        @type arg: L{escript.Data}
3707        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3708                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3709                 the eigenvector coresponding to the i-th eigenvalue.
3710        @rtype: L{tuple} of L{escript.Data}.
3711        @note: The dimension is restricted to 3.
3712        """
3713        if isinstance(arg,numarray.NumArray):
3714          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3715        elif isinstance(arg,escript.Data):
3716          return arg._eigenvalues_and_eigenvectors()
3717        elif isinstance(arg,Symbol):
3718          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3719        elif isinstance(arg,float):
3720          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3721        elif isinstance(arg,int):
3722          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3723        else:
3724          raise TypeError,"eigenvalues: Unknown argument type."
3725  #=======================================================  #=======================================================
3726  #  Binary operations:  #  Binary operations:
3727  #=======================================================  #=======================================================
# Line 3056  class Add_Symbol(DependendSymbol): Line 3765  class Add_Symbol(DependendSymbol):
3765         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3766         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3767         """         """
3768         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3769         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3770         if not sh0==sh1:         if not sh0==sh1:
3771            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3772         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 3072  class Add_Symbol(DependendSymbol): Line 3781  class Add_Symbol(DependendSymbol):
3781        @type format: C{str}        @type format: C{str}
3782        @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.
3783        @rtype: C{str}        @rtype: C{str}
3784        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3785        """        """
3786        if format=="str" or format=="text":        if format=="str" or format=="text":
3787           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 3100  class Add_Symbol(DependendSymbol): Line 3809  class Add_Symbol(DependendSymbol):
3809              raise TypeError,"%s: new value is not appropriate."%str(self)              raise TypeError,"%s: new value is not appropriate."%str(self)
3810        else:        else:
3811           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
3812           return add(args[0],args[1])           out=add(args[0],args[1])
3813             return out
3814    
3815     def diff(self,arg):     def diff(self,arg):
3816        """        """
# Line 3131  def mult(arg0,arg1): Line 3841  def mult(arg0,arg1):
3841         """         """
3842         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3843         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3844            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3845         else:         else:
3846            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3847                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3155  class Mult_Symbol(DependendSymbol): Line 3865  class Mult_Symbol(DependendSymbol):
3865         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3866         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3867         """         """
3868         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3869         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3870         if not sh0==sh1:         if not sh0==sh1:
3871            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3872         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 3171  class Mult_Symbol(DependendSymbol): Line 3881  class Mult_Symbol(DependendSymbol):
3881        @type format: C{str}        @type format: C{str}
3882        @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.
3883        @rtype: C{str}        @rtype: C{str}
3884        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3885        """        """
3886        if format=="str" or format=="text":        if format=="str" or format=="text":
3887           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3231  def quotient(arg0,arg1): Line 3941  def quotient(arg0,arg1):
3941         """         """
3942         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3943         if testForZero(args[0]):         if testForZero(args[0]):
3944            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3945         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3946            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3947               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3260  class Quotient_Symbol(DependendSymbol): Line 3970  class Quotient_Symbol(DependendSymbol):
3970         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3971         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3972         """         """
3973         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3974         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3975         if not sh0==sh1:         if not sh0==sh1:
3976            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3977         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 3276  class Quotient_Symbol(DependendSymbol): Line 3986  class Quotient_Symbol(DependendSymbol):
3986        @type format: C{str}        @type format: C{str}
3987        @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.
3988        @rtype: C{str}        @rtype: C{str}
3989        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3990        """        """
3991        if format=="str" or format=="text":        if format=="str" or format=="text":
3992           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3337  def power(arg0,arg1): Line 4047  def power(arg0,arg1):
4047         """         """
4048         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4049         if testForZero(args[0]):         if testForZero(args[0]):
4050            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4051         elif testForZero(args[1]):         elif testForZero(args[1]):
4052            return numarray.ones(args[0],numarray.Float)            return numarray.ones(getShape(args[1]),numarray.Float64)
4053         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4054            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4055         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 3362  class Power_Symbol(DependendSymbol): Line 4072  class Power_Symbol(DependendSymbol):
4072         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4073         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4074         """         """
4075         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4076         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4077         if not sh0==sh1:         if not sh0==sh1:
4078            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4079         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3380  class Power_Symbol(DependendSymbol): Line 4090  class Power_Symbol(DependendSymbol):
4090        @type format: C{str}        @type format: C{str}
4091        @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.
4092        @rtype: C{str}        @rtype: C{str}
4093        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4094        """        """
4095        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4096           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3462  def minimum(*args): Line 4172  def minimum(*args):
4172            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4173      return out      return out
4174    
4175  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4176      """      """
4177      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4178    
4179      @param arg: argument      @param arg: argument
4180      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4181      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4182      @type arg: C{float}      @type minval: C{float} or C{None}
4183      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4184      @type arg: C{float}      @type maxval: C{float} or C{None}
4185      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4186               less then maxval are unchanged.               less then maxval are unchanged.
4187      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4188      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4189      """      """
4190      if minval>maxval:      if not minval==None and not maxval==None:
4191         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4192      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4193        if minval == None:
4194            tmp=arg
4195        else:
4196            tmp=maximum(minval,arg)
4197        if maxval == None:
4198            return tmp
4199        else:
4200            return minimum(tmp,maxval)
4201    
4202        
4203  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3496  def inner(arg0,arg1): Line 4214  def inner(arg0,arg1):
4214      @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}
4215      @param arg1: second argument      @param arg1: second argument
4216      @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}
4217      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4218      @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
4219      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4220      """      """
4221      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4222      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4223      if not sh0==sh1:      if not sh0==sh1:
4224          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4225      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4226    
4227    def outer(arg0,arg1):
4228        """
4229        the outer product of the two argument:
4230    
4231        out[t,s]=arg0[t]*arg1[s]
4232    
4233        where
4234    
4235            - s runs through arg0.Shape
4236            - t runs through arg1.Shape
4237    
4238        @param arg0: first argument
4239        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4240        @param arg1: second argument
4241        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4242        @return: the outer product of arg0 and arg1 at each data point
4243        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4244        """
4245        return generalTensorProduct(arg0,arg1,axis_offset=0)
4246    
4247  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4248      """      """
4249        see L{matrix_mult}
4250        """
4251        return matrix_mult(arg0,arg1)
4252    
4253    def matrix_mult(arg0,arg1):
4254        """
4255      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4256    
4257      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4258    
4259            or      or
4260    
4261      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4262    
4263      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.
4264    
4265      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4266      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3526  def matrixmult(arg0,arg1): Line 4270  def matrixmult(arg0,arg1):
4270      @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
4271      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4272      """      """
4273      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4274      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4275      if not len(sh0)==2 :      if not len(sh0)==2 :
4276          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4277      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4278          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4279      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4280    
4281  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4282      """      """
4283      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  
4284      """      """
4285      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
   
4286    
4287  def tensormult(arg0,arg1):  def tensor_mult(arg0,arg1):
4288      """      """
4289      the tensor product of the two argument:      the tensor product of the two argument:
   
4290            
4291      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4292    
4293      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4294    
4295                   or      or
4296    
4297      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4298    
# Line 3571  def tensormult(arg0,arg1): Line 4301  def tensormult(arg0,arg1):
4301    
4302      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]
4303                                
4304                   or      or
4305    
4306      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]
4307    
4308                   or      or
4309    
4310      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]
4311    
4312      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  
4313      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.
4314    
4315      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4316      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3589  def tensormult(arg0,arg1): Line 4319  def tensormult(arg0,arg1):
4319      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4320      @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
4321      """      """
4322      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4323      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4324      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4325         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4326      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):
4327         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4328      else:      else:
4329          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4330    
4331  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4332      """      """
# Line 3604  def generalTensorProduct(arg0,arg1,axis_ Line 4334  def generalTensorProduct(arg0,arg1,axis_
4334    
4335      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4336    
4337      where s runs through arg0.Shape[:arg0.Rank-axis_offset]      where
           r runs trough arg0.Shape[:axis_offset]  
           t runs through arg1.Shape[axis_offset:]  
4338    
4339      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]
4340      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4341            - t runs through arg1.Shape[axis_offset:]
4342    
4343      @param arg0: first argument      @param arg0: first argument
4344      @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}
4345      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4346      @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}
4347      @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.
4348      @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 3626  def generalTensorProduct(arg0,arg1,axis_ Line 4355  def generalTensorProduct(arg0,arg1,axis_
4355             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4356         else:         else:
4357             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4358                 raise ValueError,"generalTensorProduct: dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_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)
4359             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4360             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4361             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 3636  def generalTensorProduct(arg0,arg1,axis_ Line 4365  def generalTensorProduct(arg0,arg1,axis_
4365             for i in sh1[:axis_offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4366             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4367             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4368             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4369             for i0 in range(d0):             for i0 in range(d0):
4370                      for i1 in range(d1):                      for i1 in range(d1):
4371                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
# Line 3652  def generalTensorProduct(arg0,arg1,axis_ Line 4381  def generalTensorProduct(arg0,arg1,axis_
4381                                    
4382  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4383     """     """
4384     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4385     """     """
4386     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4387         """         """
4388         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4389    
4390         @param arg0: numerator         @param arg0: first argument
4391         @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}.
4392         @param arg1: denominator         @param arg1: second argument
4393         @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}.
4394         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4395         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4396         """         """
4397         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4398         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4399         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4400         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4401         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 3685  class GeneralTensorProduct_Symbol(Depend Line 4414  class GeneralTensorProduct_Symbol(Depend
4414        @type format: C{str}        @type format: C{str}
4415        @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.
4416        @rtype: C{str}        @rtype: C{str}
4417        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4418        """        """
4419        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4420           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
# Line 3713  class GeneralTensorProduct_Symbol(Depend Line 4442  class GeneralTensorProduct_Symbol(Depend
4442           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4443           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4444    
4445  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4446      "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!!!"
4447      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4448      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
4449      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  def transposed_matrix_mult(arg0,arg1):
4450      shape10=arg1.getShape()[:axis_offset]      """
4451      shape1=arg1.getShape()[axis_offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4452      if not shape01==shape10:  
4453          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)      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
4454    
4455      # whatr function space should be used? (this here is not good!)      or
4456      fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
4457      # create return value:      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4458      out=escript.Data(0.,tuple(shape0+shape1),fs)  
4459      #      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4460      s0=[[]]  
4461      for k in shape0:      The first dimension of arg0 and arg1 must match.
4462            s=[]  
4463            for j in s0:      @param arg0: first argument of rank 2
4464                  for i in range(k): s.append(j+[slice(i,i)])      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4465            s0=s      @param arg1: second argument of at least rank 1
4466      s1=[[]]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4467      for k in shape1:      @return: the product of the transposed of arg0 and arg1 at each data point
4468            s=[]      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4469            for j in s1:      @raise ValueError: if the shapes of the arguments are not appropriate
4470                  for i in range(k): s.append(j+[slice(i,i)])      """
4471            s1=s      sh0=getShape(arg0)
4472      s01=[[]]      sh1=getShape(arg1)
4473      for k in shape01:      if not len(sh0)==2 :
4474            s=[]          raise ValueError,"first argument must have rank 2"
4475            for j in s01:      if not len(sh1)==2 and not len(sh1)==1:
4476                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"second argument must have rank 1 or 2"
4477            s01=s      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4478    
4479      for i0 in s0:  def transposed_tensor_mult(arg0,arg1):
4480         for i1 in s1:      """
4481           s=escript.Scalar(0.,fs)      the tensor product of the transposed of the first and the second argument
4482           for i01 in s01:      
4483              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      for arg0 of rank 2 this is
4484           out.__setitem__(tuple(i0+i1),s)  
4485      return out      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4486    
4487        or
4488    
4489        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4490    
4491      
4492        and for arg0 of rank 4 this is
4493    
4494        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4495                  
4496        or
4497    
4498        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4499    
4500        or
4501    
4502        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4503    
4504        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4505        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4506    
4507        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4508    
4509        @param arg0: first argument of rank 2 or 4
4510        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4511        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4512        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4513        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4514        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4515        """
4516        sh0=getShape(arg0)
4517        sh1=getShape(arg1)
4518        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4519           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4520        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4521           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4522        else:
4523            raise ValueError,"first argument must have rank 2 or 4"
4524    
4525    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4526        """
4527        generalized tensor product of transposed of arg0 and arg1:
4528    
4529        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4530    
4531        where
4532    
4533            - s runs through arg0.Shape[axis_offset:]
4534            - r runs trough arg0.Shape[:axis_offset]
4535            - t runs through arg1.Shape[axis_offset:]
4536    
4537        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4538        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4539    
4540        @param arg0: first argument
4541        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4542        @param arg1: second argument
4543        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4544        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4545        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4546        """
4547        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4548        arg0,arg1=matchType(arg0,arg1)
4549        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4550        if isinstance(arg0,numarray.NumArray):
4551           if isinstance(arg1,Symbol):
4552               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4553           else:
4554               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4555                   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)
4556               arg0_c=arg0.copy()
4557               arg1_c=arg1.copy()
4558               sh0,sh1=arg0.shape,arg1.shape
4559               d0,d1,d01=1,1,1
4560               for i in sh0[axis_offset:]: d0*=i
4561               for i in sh1[axis_offset:]: d1*=i
4562               for i in sh0[:axis_offset]: d01*=i
4563               arg0_c.resize((d01,d0))
4564               arg1_c.resize((d01,d1))
4565               out=numarray.zeros((d0,d1),numarray.Float64)
4566               for i0 in range(d0):
4567                        for i1 in range(d1):
4568                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4569               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4570               return out
4571        elif isinstance(arg0,escript.Data):
4572           if isinstance(arg1,Symbol):
4573               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4574           else:
4575               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4576        else:      
4577           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4578                    
4579    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4580       """
4581       Symbol representing the general tensor product of the transposed of arg0 and arg1
4582       """
4583       def __init__(self,arg0,arg1,axis_offset=0):
4584           """
4585           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4586    
4587           @param arg0: first argument
4588           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4589           @param arg1: second argument
4590           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4591           @raise ValueError: inconsistent dimensions of arguments.
4592           @note: if both arguments have a spatial dimension, they must equal.
4593           """
4594           sh_arg0=getShape(arg0)
4595           sh_arg1=getShape(arg1)
4596           sh01=sh_arg0[:axis_offset]
4597           sh10=sh_arg1[:axis_offset]
4598           sh0=sh_arg0[axis_offset:]
4599           sh1=sh_arg1[axis_offset:]
4600           if not sh01==sh10:
4601               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)
4602           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4603    
4604       def getMyCode(self,argstrs,format="escript"):
4605          """
4606          returns a program code that can be used to evaluate the symbol.
4607    
4608          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4609          @type argstrs: C{list} of length 2 of C{str}.
4610          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4611          @type format: C{str}
4612          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4613          @rtype: C{str}
4614          @raise NotImplementedError: if the requested format is not available
4615          """
4616          if format=="escript" or format=="str" or format=="text":
4617             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4618          else:
4619             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4620    
4621       def substitute(self,argvals):
4622          """
4623          assigns new values to symbols in the definition of the symbol.
4624          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4625    
4626          @param argvals: new values assigned to symbols
4627          @type argvals: C{dict} with keywords of type L{Symbol}.
4628          @return: result of the substitution process. Operations are executed as much as possible.
4629          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4630          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4631          """
4632          if argvals.has_key(self):
4633             arg=argvals[self]
4634             if self.isAppropriateValue(arg):
4635                return arg
4636             else:
4637                raise TypeError,"%s: new value is not appropriate."%str(self)
4638          else:
4639             args=self.getSubstitutedArguments(argvals)
4640             return generalTransposedTensorProduct(args[0],args[1],args[2])
4641    
4642    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4643        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4644        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4645    
4646    def matrix_transposed_mult(arg0,arg1):
4647        """
4648        matrix-transposed(matrix) product of the two argument:
4649    
4650        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4651    
4652        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4653    
4654        The last dimensions of arg0 and arg1 must match.
4655    
4656        @param arg0: first argument of rank 2
4657        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4658        @param arg1: second argument of rank 2
4659        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4660        @return: the product of arg0 and the transposed of arg1 at each data point
4661        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4662        @raise ValueError: if the shapes of the arguments are not appropriate
4663        """
4664        sh0=getShape(arg0)
4665        sh1=getShape(arg1)
4666        if not len(sh0)==2 :
4667            raise ValueError,"first argument must have rank 2"
4668        if not len(sh1)==2 and not len(sh1)==1:
4669            raise ValueError,"second argument must have rank 1 or 2"
4670        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4671    
4672    def tensor_transposed_mult(arg0,arg1):
4673        """
4674        the tensor product of the first and the transpose of the second argument
4675        
4676        for arg0 of rank 2 this is
4677    
4678        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4679    
4680        and for arg0 of rank 4 this is
4681    
4682        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4683                  
4684        or
4685    
4686        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4687    
4688        In the first case the the second dimension of arg0 and arg1 must match and  
4689        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4690    
4691        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4692    
4693        @param arg0: first argument of rank 2 or 4
4694        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4695        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4696        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4697        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4698        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4699        """
4700        sh0=getShape(arg0)
4701        sh1=getShape(arg1)
4702        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4703           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4704        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4705           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4706        else:
4707            raise ValueError,"first argument must have rank 2 or 4"
4708    
4709    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4710        """
4711        generalized tensor product of transposed of arg0 and arg1:
4712    
4713        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4714    
4715        where
4716    
4717            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4718            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4719            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4720    
4721        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4722        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4723    
4724        @param arg0: first argument
4725        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4726        @param arg1: second argument
4727        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4728        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4729        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4730        """
4731        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4732        arg0,arg1=matchType(arg0,arg1)
4733        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4734        if isinstance(arg0,numarray.NumArray):
4735           if isinstance(arg1,Symbol):
4736               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4737           else:
4738               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4739                   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)
4740               arg0_c=arg0.copy()
4741               arg1_c=arg1.copy()
4742               sh0,sh1=arg0.shape,arg1.shape
4743               d0,d1,d01=1,1,1
4744               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4745               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4746               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4747               arg0_c.resize((d0,d01))
4748               arg1_c.resize((d1,d01))
4749               out=numarray.zeros((d0,d1),numarray.Float64)
4750               for i0 in range(d0):
4751                        for i1 in range(d1):
4752                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4753               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4754               return out
4755        elif isinstance(arg0,escript.Data):
4756           if isinstance(arg1,Symbol):
4757               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4758           else:
4759               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4760        else:      
4761           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4762                    
4763    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4764       """
4765       Symbol representing the general tensor product of arg0 and the transpose of arg1
4766       """
4767       def __init__(self,arg0,arg1,axis_offset=0):
4768           """
4769           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4770    
4771           @param arg0: first argument
4772           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4773           @param arg1: second argument
4774           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4775           @raise ValueError: inconsistent dimensions of arguments.
4776           @note: if both arguments have a spatial dimension, they must equal.
4777           """
4778           sh_arg0=getShape(arg0)
4779           sh_arg1=getShape(arg1)
4780           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4781           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4782           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4783           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4784           if not sh01==sh10:
4785               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)
4786           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4787    
4788       def getMyCode(self,argstrs,format="escript"):
4789          """
4790          returns a program code that can be used to evaluate the symbol.
4791    
4792          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4793          @type argstrs: C{list} of length 2 of C{str}.
4794          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4795          @type format: C{str}
4796          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4797          @rtype: C{str}
4798          @raise NotImplementedError: if the requested format is not available
4799          """
4800          if format=="escript" or format=="str" or format=="text":
4801             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4802          else:
4803             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4804    
4805       def substitute(self,argvals):
4806          """
4807          assigns new values to symbols in the definition of the symbol.
4808          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4809    
4810          @param argvals: new values assigned to symbols
4811          @type argvals: C{dict} with keywords of type L{Symbol}.
4812          @return: result of the substitution process. Operations are executed as much as possible.
4813          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4814          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4815          """
4816          if argvals.has_key(self):
4817             arg=argvals[self]
4818             if self.isAppropriateValue(arg):
4819                return arg
4820             else:
4821                raise TypeError,"%s: new value is not appropriate."%str(self)
4822          else:
4823             args=self.getSubstitutedArguments(argvals)
4824             return generalTensorTransposedProduct(args[0],args[1],args[2])
4825    
4826    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4827        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4828        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4829    
4830  #=========================================================  #=========================================================
4831  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
# Line 3776  def grad(arg,where=None): Line 4847  def grad(arg,where=None):
4847                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4848      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4849      @return: gradient of arg.      @return: gradient of arg.
4850      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4851      """      """
4852      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4853         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3804  class Grad_Symbol(DependendSymbol): Line 4875  class Grad_Symbol(DependendSymbol):
4875        d=arg.getDim()        d=arg.getDim()
4876        if d==None:        if d==None:
4877           raise ValueError,"argument must have a spatial dimension"           raise ValueError,"argument must have a spatial dimension"
4878        super(Grad_Symbol,self).__init__(args=[arg,where],shape=tuple(list(arg.getShape()).extend(d)),dim=d)        super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4879    
4880     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4881        """        """
# Line 3816  class Grad_Symbol(DependendSymbol): Line 4887  class Grad_Symbol(DependendSymbol):
4887        @type format: C{str}        @type format: C{str}
4888        @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.
4889        @rtype: C{str}        @rtype: C{str}
4890        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4891        """        """
4892        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4893           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 3869  def integrate(arg,where=None): Line 4940  def integrate(arg,where=None):
4940                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4941      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4942      @return: integral of arg.      @return: integral of arg.
4943      @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}      @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4944      """      """
4945      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4946         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
# Line 3880  def integrate(arg,where=None): Line 4951  def integrate(arg,where=None):
4951         else:         else:
4952            return arg._integrate()            return arg._integrate()
4953      else:      else:
4954        raise TypeError,"integrate: Unknown argument type."         arg2=escript.Data(arg,where)
4955           if arg2.getRank()==0:
4956              return arg2._integrate()[0]
4957           else:
4958              return arg2._integrate()
4959    
4960  class Integrate_Symbol(DependendSymbol):  class Integrate_Symbol(DependendSymbol):
4961     """     """
# Line 3907  class Integrate_Symbol(DependendSymbol): Line 4982  class Integrate_Symbol(DependendSymbol):
4982        @type format: C{str}        @type format: C{str}
4983        @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.
4984        @rtype: C{str}        @rtype: C{str}
4985        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4986        """        """
4987        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4988           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 3952  class Integrate_Symbol(DependendSymbol): Line 5027  class Integrate_Symbol(DependendSymbol):
5027    
5028  def interpolate(arg,where):  def interpolate(arg,where):
5029      """      """
5030      interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where. If the argument C{arg} has the requested function space
5031        C{where} no interpolation is performed and C{arg} is returned.
5032    
5033      @param arg: interpolant      @param arg: interpolant
5034      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
5035      @param where: FunctionSpace to be interpolated to      @param where: FunctionSpace to be interpolated to
5036      @type where: L{escript.FunctionSpace}      @type where: L{escript.FunctionSpace}
5037      @return: interpolated argument      @return: interpolated argument
5038      @rtype:  C{escript.Data} or L{Symbol}      @rtype: C{escript.Data} or L{Symbol}
5039      """      """
5040      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5041         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
5042        elif isinstance(arg,escript.Data):
5043           if arg.isEmpty():
5044              return arg
5045           elif where == arg.getFunctionSpace():
5046              return arg
5047           else:
5048              return escript.Data(arg,where)
5049      else:      else:
5050         return escript.Data(arg,where)         return escript.Data(arg,where)
5051    
# Line 3990  class Interpolate_Symbol(DependendSymbol Line 5073  class Interpolate_Symbol(DependendSymbol
5073        @type format: C{str}        @type format: C{str}
5074        @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.
5075        @rtype: C{str}        @rtype: C{str}
5076        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
5077        """        """
5078        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
5079           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4043  def div(arg,where=None): Line 5126  def div(arg,where=None):
5126                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
5127      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
5128      @return: divergence of arg.      @return: divergence of arg.
5129      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5130      """      """
5131      if not arg.getShape()==(arg.getDim(),):      if isinstance(arg,Symbol):
5132        raise ValueError,"div: expected shape is (%s,)"%arg.getDim()          dim=arg.getDim()
5133        elif isinstance(arg,escript.Data):
5134            dim=arg.getDomain().getDim()
5135        else:
5136            raise TypeError,"div: argument type not supported"
5137        if not arg.getShape()==(dim,):
5138          raise ValueError,"div: expected shape is (%s,)"%dim
5139      return trace(grad(arg,where))      return trace(grad(arg,where))
5140    
5141  def jump(arg,domain=None):  def jump(arg,domain=None):
# Line 4059  def jump(arg,domain=None): Line 5148  def jump(arg,domain=None):
5148                     the domain of arg is used. If arg is a L{Symbol} the domain must be present.                     the domain of arg is used. If arg is a L{Symbol} the domain must be present.
5149      @type domain: C{None} or L{escript.Domain}      @type domain: C{None} or L{escript.Domain}
5150      @return: jump of arg      @return: jump of arg
5151      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5152      """      """
5153      if domain==None: domain=arg.getDomain()      if domain==None: domain=arg.getDomain()
5154      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
 #=============================  
 #  
 # wrapper for various functions: if the argument has attribute the function name  
 # as an argument it calls the corresponding methods. Otherwise the corresponding  
 # numarray function is called.  
5155    
5156  # functions involving the underlying Domain:  def L2(arg):
5157        """
5158        returns the L2 norm of arg at where
5159        
5160        @param arg: function which L2 to be calculated.
5161        @type arg: L{escript.Data} or L{Symbol}
5162        @return: L2 norm of arg.
5163        @rtype: L{float} or L{Symbol}
5164        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5165        """
5166        return sqrt(integrate(inner(arg,arg)))
5167    
5168    def getClosestValue(arg,origin=0):
5169        """
5170        returns the value in arg which is closest to origin
5171        
5172        @param arg: function
5173        @type arg: L{escript.Data} or L{Symbol}
5174        @param origin: reference value
5175        @type origin: C{float} or L{escript.Data}
5176        @return: value in arg closest to origin
5177        @rtype: L{numarray.NumArray} or L{Symbol}
5178        """
5179        return arg.getValueOfGlobalDataPoint(*(length(arg-origin).minGlobalDataPoint()))
5180    
5181  def transpose(arg,axis=None):  def normalize(arg,zerolength=0):
5182        """
5183        returns normalized version of arg (=arg/length(arg))
5184        
5185        @param arg: function
5186        @type arg: L{escript.Data} or L{Symbol}
5187        @param zerolength: realitive tolerance for arg == 0.
5188        @type zerolength: C{float}
5189        @return: normalized arg where arg is non zero and zero elsewhere
5190        @rtype: L{escript.Data} or L{Symbol}
5191      """      """
5192      Returns the transpose of the Data object arg.      l=length(arg)
5193        m=whereZero(l,zerolength*Lsup(l))
5194        mm=1-m
5195        return arg*(mm/(l*mm+m))
5196    
5197      @param arg:  def deviatoric(arg):
5198      """      """
5199      if axis==None:      returns the deviatoric version of arg
5200         r=0      """
5201         if hasattr(arg,"getRank"): r=arg.getRank()      return arg-(trace(arg)/trace(kronecker(arg.getDomain())))*kronecker(arg.getDomain())
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
     if isinstance(arg,Symbol):  
        return Transpose_Symbol(arg,axis=r)  
     if isinstance(arg,escript.Data):  
        # hack for transpose  
        r=arg.getRank()  
        if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"  
        s=arg.getShape()  
        out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())  
        for i in range(s[0]):  
           for j in range(s[1]):  
              out[j,i]=arg[i,j]  
        return out  
        # end hack for transpose  
        return arg.transpose(axis)  
     else:  
        return numarray.transpose(arg,axis=axis)  
5202    
5203    def diameter(domain):
5204        """
5205        returns the diameter of a domain
5206    
5207        @param domain: a domain
5208        @type domain: L{escript.Domain}
5209        @rtype: C{float}
5210        """
5211        x=domain.getX()
5212        out=0.
5213        for i in xrange(domain.getDim()):
5214           x_i=x[i]
5215           out=max(out,sup(x_i)-inf(x_i))
5216        return out
5217    #=============================
5218    #
5219    
5220  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5221      """      """
5222      resorts the component of arg according to index      resorts the component of arg according to index
5223    
5224      """      """
5225      pass      raise NotImplementedError
 #  
 # $Log: util.py,v $  
 # Revision 1.14.2.16  2005/10/19 06:09:57  gross  
 # new versions of saveVTK and saveDX which allow to write several Data objects into a single file  
 #  
 # Revision 1.14.2.15  2005/10/19 03:25:27  jgs  
 # numarray uses log10 for log, and log for ln - log()/ln() updated to reflect this  
 #  
 # Revision 1.14.2.14  2005/10/19 02:10:23  jgs  
 # fixed ln() to call Data::ln  
 #  
 # Revision 1.14.2.13  2005/09/12 03:32:14  gross  
 # test_visualiztion has been aded to mk  
 #  
 # Revision 1.14.2.12  2005/09/09 01:56:24  jgs  
 # added implementations of acos asin atan sinh cosh tanh asinh acosh atanh  
 # and some associated testing  
 #  
 # Revision 1.14.2.11  2005/09/08 08:28:39  gross  
 # some cleanup in savevtk  
 #  
 # Revision 1.14.2.10  2005/09/08 00:25:32  gross  
 # test for finley mesh generators added  
 #  
 # Revision 1.14.2.9  2005/09/07 10:32:05  gross  
 # Symbols removed from util and put into symmbols.py.  
 #  
 # Revision 1.14.2.8  2005/08/26 05:06:37  cochrane  
 # Corrected errors in docstrings.  Improved output formatting of docstrings.  
 # Other minor improvements to code and docs (eg spelling etc).  
 #  
 # Revision 1.14.2.7  2005/08/26 04:45:40  cochrane  
 # Fixed and tidied markup and docstrings.  Some *Symbol classes were defined  
 # as functions, so changed them to classes (hopefully this was the right thing  
 # to do).  
 #  
 # Revision 1.14.2.6  2005/08/26 04:30:13  gross  
 # gneric unit testing for linearPDE  
 #  
 # Revision 1.14.2.5  2005/08/24 02:02:52  gross  
 # jump function added  
 #  
 # Revision 1.14.2.4  2005/08/18 04:39:32  gross  
 # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now  
 #  
 # Revision 1.14.2.3  2005/08/03 09:55:33  gross  
 # ContactTest is passing now./mk install!  
 #  
 # Revision 1.14.2.2  2005/08/02 03:15:14  gross  
 # bug inb trace fixed!  
 #  
 # Revision 1.14.2.1  2005/07/29 07:10:28  gross  
 # new functions in util and a new pde type in linearPDEs  
 #  
 # Revision 1.2.2.21  2005/07/28 04:19:23  gross  
 # new functions maximum and minimum introduced.  
 #  
 # Revision 1.2.2.20  2005/07/25 01:26:27  gross  
 # bug in inner fixed  
 #  
 # Revision 1.2.2.19  2005/07/21 04:01:28  jgs  
 # minor comment fixes  
 #  
 # Revision 1.2.2.18  2005/07/21 01:02:43  jgs  
 # commit ln() updates to development branch version  
 #  
 # Revision 1.12  2005/07/20 06:14:58  jgs  
 # added ln(data) style wrapper for data.ln() - also added corresponding  
 # implementation of Ln_Symbol class (not sure if this is right though)  
 #  
 # Revision 1.11  2005/07/08 04:07:35  jgs  
 # Merge of development branch back to main trunk on 2005-07-08  
 #  
 # Revision 1.10  2005/06/09 05:37:59  jgs  
 # Merge of development branch back to main trunk on 2005-06-09  
 #  
 # Revision 1.2.2.17  2005/07/07 07:28:58  gross  
 # some stuff added to util.py to improve functionality  
 #  
 # Revision 1.2.2.16  2005/06/30 01:53:55  gross  
 # a bug in coloring fixed  
 #  
 # Revision 1.2.2.15  2005/06/29 02:36:43  gross  
 # Symbols have been introduced and some function clarified. needs much more work  
 #  
 # Revision 1.2.2.14  2005/05/20 04:05:23  gross  
 # some work on a darcy flow started  
 #  
 # Revision 1.2.2.13  2005/03/16 05:17:58  matt  
 # Implemented unit(idx, dim) to create cartesian unit basis vectors to  
 # complement kronecker(dim) function.  
 #  
 # Revision 1.2.2.12  2005/03/10 08:14:37  matt  
 # Added non-member Linf utility function to complement Data::Linf().  
 #  
 # Revision 1.2.2.11  2005/02/17 05:53:25  gross  
 # some bug in saveDX fixed: in fact the bug was in  
 # DataC/getDataPointShape  
 #  
 # Revision 1.2.2.10  2005/01/11 04:59:36  gross  
 # automatic interpolation in integrate switched off  
 #  
 # Revision 1.2.2.9  2005/01/11 03:38:13  gross  
 # Bug in Data.integrate() fixed for the case of rank 0. The problem is not totallly resolved as the method should return a scalar rather than a numarray object in the case of rank 0. This problem is fixed by the util.integrate wrapper.  
 #  
 # Revision 1.2.2.8  2005/01/05 04:21:41  gross  
 # FunctionSpace checking/matchig in slicing added  
 #  
 # Revision 1.2.2.7  2004/12/29 05:29:59  gross  
 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()  
 #  
 # Revision 1.2.2.6  2004/12/24 06:05:41  gross  
 # some changes in linearPDEs to add AdevectivePDE  
 #  
 # Revision 1.2.2.5  2004/12/17 00:06:53  gross  
 # mk sets ESYS_ROOT is undefined  
 #  
 # Revision 1.2.2.4  2004/12/07 03:19:51  gross  
 # options for GMRES and PRES20 added  
 #  
 # Revision 1.2.2.3  2004/12/06 04:55:18  gross  
 # function wraper extended  
 #  
 # Revision 1.2.2.2  2004/11/22 05:44:07  gross  
 # a few more unitary functions have been added but not implemented in Data yet  
 #  
 # Revision 1.2.2.1  2004/11/12 06:58:15  gross  
 # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry  
 #  
 # Revision 1.2  2004/10/27 00:23:36  jgs  
 # fixed minor syntax error  
 #  
 # Revision 1.1.1.1  2004/10/26 06:53:56  jgs  
 # initial import of project esys2  
 #  
 # Revision 1.1.2.3  2004/10/26 06:43:48  jgs  
 # committing Lutz's and Paul's changes to brach jgs  
 #  
 # Revision 1.1.4.1  2004/10/20 05:32:51  cochrane  
 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.  
 #  
 # Revision 1.1  2004/08/05 03:58:27  gross  
 # Bug in Assemble_NodeCoordinates fixed  
 #  
 #  
5226    
5227  # vim: expandtab shiftwidth=4:  def showEscriptParams():
5228        """
5229        Display the parameters escript recognises with an explanation.
5230        """
5231        p=listEscriptParams()
5232        for name,desc in p:
5233        print name+':\t'+desc

Legend:
Removed from v.429  
changed lines
  Added in v.2142

  ViewVC Help
Powered by ViewVC 1.1.26