/[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 587 by gross, Fri Mar 10 02:26:50 2006 UTC revision 1357 by gross, Thu Nov 29 08:39:32 2007 UTC
# Line 1  Line 1 
1    #
2  # $Id$  # $Id$
3  #  #
4  #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  #######################################################
5    #
6    #           Copyright 2003-2007 by ACceSS MNRF
7    #       Copyright 2007 by University of Queensland
8  #  #
9  #   This software is the property of ACcESS.  No part of this code  #                http://esscc.uq.edu.au
10  #   may be copied in any form or by any means without the expressed written  #        Primary Business: Queensland, Australia
11  #   consent of ACcESS.  Copying, use or modification of this software  #  Licensed under the Open Software License version 3.0
12  #   by any unauthorised person is illegal unless that  #     http://www.opensource.org/licenses/osl-3.0.php
13  #   person has a software license agreement with ACcESS.  #
14    #######################################################
15  #  #
16    
17  """  """
18  Utility functions for escript  Utility functions for escript
19    
 @remark:  This module is under construction and is still tested!!!  
   
20  @var __author__: name of author  @var __author__: name of author
21  @var __licence__: licence agreement  @var __copyright__: copyrights
22    @var __license__: licence agreement
23  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
24  @var __version__: version  @var __version__: version
25  @var __date__: date of the version  @var __date__: date of the version
26    @var EPSILON: smallest positive value with 1.<1.+EPSILON
27  """  """
28                                                                                                                                                                                                                                                                                                                                                                                                            
29  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
30  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
31                        http://www.access.edu.au
32                    Primary Business: Queensland, Australia"""
33    __license__="""Licensed under the Open Software License version 3.0
34                 http://www.opensource.org/licenses/osl-3.0.php"""
35  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
36  __version__="$Revision$"  __version__="$Revision$"
37  __date__="$Date$"  __date__="$Date$"
# Line 30  __date__="$Date$" Line 39  __date__="$Date$"
39    
40  import math  import math
41  import numarray  import numarray
 import numarray.linear_algebra  
42  import escript  import escript
43  import os  import os
44    from esys.escript import C_GeneralTensorProduct
45    from esys.escript import getVersion
46    
47  # missing tests:  #=========================================================
48    #   some helpers:
49    #=========================================================
50    def getEpsilon():
51         #     ------------------------------------------------------------------
52         #     Compute EPSILON, the machine precision.  The call to daxpp is
53         #     inTENded to fool compilers that use extra-length registers.
54         #     31 Map 1999: Hardwire EPSILON so the debugger can step thru easily.
55         #     ------------------------------------------------------------------
56         eps    = 2.**(-12)
57         p=2.
58         while p>1.:
59                eps/=2.
60                p=1.+eps
61         return eps*2.
62    
63  # def pokeShape(arg):  EPSILON=getEpsilon()
 # def pokeDim(arg):  
 # def commonShape(arg0,arg1):  
 # def commonDim(*args):  
 # def testForZero(arg):  
 # def matchType(arg0=0.,arg1=0.):  
 # def matchShape(arg0,arg1):  
64    
65  # def reorderComponents(arg,index):  def getTagNames(domain):
66        """
67        returns a list of the tag names used by the domain
68    
69  #      
70  # slicing: get      @param domain: a domain object
71  #          set      @type domain: L{escript.Domain}
72  #      @return: a list of the tag name used by the domain.
73  # and derivatives      @rtype: C{list} of C{str}
74        """
75        return [n.strip() for n in domain.showTagNames().split(",") ]
76    
77    def insertTagNames(domain,**kwargs):
78        """
79        inserts tag names into the domain
80    
81        @param domain: a domain object
82        @type domain: C{escript.Domain}
83        @keyword <tag name>: tag key assigned to <tag name>
84        @type <tag name>: C{int}
85        """
86        for  k in kwargs:
87             domain.setTagMap(k,kwargs[k])
88    
89    def insertTaggedValues(target,**kwargs):
90        """
91        inserts tagged values into the tagged using tag names
92    
93        @param target: data to be filled by tagged values
94        @type target: L{escript.Data}
95        @keyword <tag name>: value to be used for <tag name>
96        @type <tag name>: C{float} or {numarray.NumArray}
97        @return: C{target}
98        @rtype: L{escript.Data}
99        """
100        for k in kwargs:
101            target.setTaggedValue(k,kwargs[k])
102        return target
103    
 #=========================================================  
 #   some helpers:  
 #=========================================================  
104  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
105      """      """
106      writes a L{Data} objects into a files using the the VTK XML file format.      writes a L{Data} objects into a files using the the VTK XML file format.
107    
108      Example:      Example::
109    
110         tmp=Scalar(..)         tmp=Scalar(..)
111         v=Vector(..)         v=Vector(..)
112         saveVTK("solution.xml",temperature=tmp,velovity=v)         saveVTK("solution.xml",temperature=tmp,velocity=v)
113    
114      tmp and v are written into "solution.xml" where tmp is named "temperature" and v is named "velovity"      tmp and v are written into "solution.xml" where tmp is named "temperature" and v is named "velocity"
115    
116      @param filename: file name of the output file      @param filename: file name of the output file
117      @type filename: C{str}      @type filename: C{str}
# Line 75  def saveVTK(filename,domain=None,**data) Line 121  def saveVTK(filename,domain=None,**data)
121      @type <name>: L{Data} object.      @type <name>: L{Data} object.
122      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
123      """      """
124      if domain==None:      new_data={}
125         for i in data.keys():      for n,d in data.items():
126            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
127                fs=d.getFunctionSpace()
128                domain2=fs.getDomain()
129                if fs == escript.Solution(domain2):
130                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
131                elif fs == escript.ReducedSolution(domain2):
132                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
133                else:
134                   new_data[n]=d
135                if domain==None: domain=domain2
136      if domain==None:      if domain==None:
137          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
138      else:      domain.saveVTK(filename,new_data)
         domain.saveVTK(filename,data)  
139    
140  def saveDX(filename,domain=None,**data):  def saveDX(filename,domain=None,**data):
141      """      """
142      writes a L{Data} objects into a files using the the DX file format.      writes a L{Data} objects into a files using the the DX file format.
143    
144      Example:      Example::
145    
146         tmp=Scalar(..)         tmp=Scalar(..)
147         v=Vector(..)         v=Vector(..)
148         saveDX("solution.dx",temperature=tmp,velovity=v)         saveDX("solution.dx",temperature=tmp,velocity=v)
149    
150      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity".      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velocity".
151    
152      @param filename: file name of the output file      @param filename: file name of the output file
153      @type filename: C{str}      @type filename: C{str}
# Line 103  def saveDX(filename,domain=None,**data): Line 157  def saveDX(filename,domain=None,**data):
157      @type <name>: L{Data} object.      @type <name>: L{Data} object.
158      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
159      """      """
160      if domain==None:      new_data={}
161         for i in data.keys():      for n,d in data.items():
162            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
163                fs=d.getFunctionSpace()
164                domain2=fs.getDomain()
165                if fs == escript.Solution(domain2):
166                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
167                elif fs == escript.ReducedSolution(domain2):
168                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
169                elif fs == escript.ContinuousFunction(domain2):
170                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
171                else:
172                   new_data[n]=d
173                if domain==None: domain=domain2
174      if domain==None:      if domain==None:
175          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
176      else:      domain.saveDX(filename,new_data)
         domain.saveDX(filename,data)  
177    
178  def kronecker(d=3):  def kronecker(d=3):
179     """     """
# Line 118  def kronecker(d=3): Line 182  def kronecker(d=3):
182     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
183     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
184     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
185     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2.
186     """     """
187     return identityTensor(d)     return identityTensor(d)
188    
# Line 154  def identityTensor(d=3): Line 218  def identityTensor(d=3):
218     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
219     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
220     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise     @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
221     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
222     """     """
223     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
224         return escript.Data(identity((d.getDim(),)),d)         return escript.Data(identity((d.getDim(),)),d)
# Line 170  def identityTensor4(d=3): Line 234  def identityTensor4(d=3):
234     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
235     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
236     @return: the object u of rank 4 with M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise     @return: the object u of rank 4 with M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise
237     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
238     """     """
239     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
240         return escript.Data(identity((d.getDim(),d.getDim())),d)         return escript.Data(identity((d.getDim(),d.getDim())),d)
# Line 188  def unitVector(i=0,d=3): Line 252  def unitVector(i=0,d=3):
252     @param d: dimension or an object that has the C{getDim} method defining the dimension     @param d: dimension or an object that has the C{getDim} method defining the dimension
253     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
254     @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise     @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise
255     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 1     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
256     """     """
257     return kronecker(d)[i]     return kronecker(d)[i]
258    
# Line 244  def inf(arg): Line 308  def inf(arg):
308    
309      @param arg: argument      @param arg: argument
310      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
311      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
312      @rtype: C{float}      @rtype: C{float}
313      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
314      """      """
# Line 263  def inf(arg): Line 327  def inf(arg):
327  #=========================================================================  #=========================================================================
328  #   some little helpers  #   some little helpers
329  #=========================================================================  #=========================================================================
330  def pokeShape(arg):  def getRank(arg):
331        """
332        identifies the rank of its argument
333    
334        @param arg: a given object
335        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
336        @return: the rank of the argument
337        @rtype: C{int}
338        @raise TypeError: if type of arg cannot be processed
339        """
340    
341        if isinstance(arg,numarray.NumArray):
342            return arg.rank
343        elif isinstance(arg,escript.Data):
344            return arg.getRank()
345        elif isinstance(arg,float):
346            return 0
347        elif isinstance(arg,int):
348            return 0
349        elif isinstance(arg,Symbol):
350            return arg.getRank()
351        else:
352          raise TypeError,"getShape: cannot identify shape"
353    def getShape(arg):
354      """      """
355      identifies the shape of its argument      identifies the shape of its argument
356    
# Line 285  def pokeShape(arg): Line 372  def pokeShape(arg):
372      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
373          return arg.getShape()          return arg.getShape()
374      else:      else:
375        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
376    
377  def pokeDim(arg):  def pokeDim(arg):
378      """      """
# Line 308  def commonShape(arg0,arg1): Line 395  def commonShape(arg0,arg1):
395      """      """
396      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.
397    
398      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
399      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
400      @return: the shape of arg0 or arg1 such that the left port equals the shape of arg0 and the right end equals the shape of arg1.      @return: the shape of arg0 or arg1 such that the left port equals the shape of arg0 and the right end equals the shape of arg1.
401      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
402      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
403      """      """
404      sh0=pokeShape(arg0)      sh0=getShape(arg0)
405      sh1=pokeShape(arg1)      sh1=getShape(arg1)
406      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
407         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
408               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
# Line 333  def commonDim(*args): Line 420  def commonDim(*args):
420      """      """
421      identifies, if possible, the spatial dimension across a set of objects which may or my not have a spatial dimension.      identifies, if possible, the spatial dimension across a set of objects which may or my not have a spatial dimension.
422    
423      @param *args: given objects      @param args: given objects
424      @return: the spatial dimension of the objects with identifiable dimension (see L{pokeDim}). If none the objects has      @return: the spatial dimension of the objects with identifiable dimension (see L{pokeDim}). If none the objects has
425               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
426      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 355  def testForZero(arg): Line 442  def testForZero(arg):
442    
443      @param arg: a given object      @param arg: a given object
444      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}
445      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
446      @rtype : C{bool}      @rtype: C{bool}
447      """      """
448      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
449         return not Lsup(arg)>0.         return not Lsup(arg)>0.
# Line 459  def matchType(arg0=0.,arg1=0.): Line 546  def matchType(arg0=0.,arg1=0.):
546    
547  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
548      """      """
549            return representations of arg0 amd arg1 which ahve the same shape
   
     If shape is not given the shape "largest" shape of args is used.  
550    
551      @param args: a given ob      @param arg0: a given object
552      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
553      @return: True if the argument is identical to zero.      @param arg1: a given object
554      @rtype: C{list} of C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
555        @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
556        @rtype: C{tuple}
557      """      """
558      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
559      sh0=pokeShape(arg0)      sh0=getShape(arg0)
560      sh1=pokeShape(arg1)      sh1=getShape(arg1)
561      if len(sh0)<len(sh):      if len(sh0)<len(sh):
562         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
563      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
# Line 494  class Symbol(object): Line 581  class Symbol(object):
581         Creates an instance of a symbol of a given shape. The symbol may depending on a list of arguments args which may be         Creates an instance of a symbol of a given shape. The symbol may depending on a list of arguments args which may be
582         symbols or any other object.         symbols or any other object.
583    
584         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
585         @type arg: C{list}         @type args: C{list}
586         @param shape: the shape         @param shape: the shape
587         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
588         @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
# Line 538  class Symbol(object): Line 625  class Symbol(object):
625         """         """
626         the shape of the symbol.         the shape of the symbol.
627    
628         @return : the shape of the symbol.         @return: the shape of the symbol.
629         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
630         """         """
631         return self.__shape         return self.__shape
# Line 547  class Symbol(object): Line 634  class Symbol(object):
634         """         """
635         the spatial dimension         the spatial dimension
636    
637         @return : the spatial dimension         @return: the spatial dimension
638         @rtype: C{int} if the dimension is defined. Otherwise C{None} is returned.         @rtype: C{int} if the dimension is defined. Otherwise C{None} is returned.
639         """         """
640         return self.__dim         return self.__dim
# Line 571  class Symbol(object): Line 658  class Symbol(object):
658         """         """
659         substitutes symbols in the arguments of this object and returns the result as a list.         substitutes symbols in the arguments of this object and returns the result as a list.
660    
661         @param argvals: L{Symbols} and their substitutes. The L{Symbol} u in the expression defining this object is replaced by argvals[u].         @param argvals: L{Symbol} and their substitutes. The L{Symbol} u in the expression defining this object is replaced by argvals[u].
662         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
663         @rtype: C{list} of objects         @rtype: C{list} of objects
664         @return: list of the object assigned to the arguments through substitution or for the arguments which are not L{Symbols} the value assigned to the argument at instantiation.         @return: list of the object assigned to the arguments through substitution or for the arguments which are not L{Symbol} the value assigned to the argument at instantiation.
665         """         """
666         out=[]         out=[]
667         for a in self.getArgument():         for a in self.getArgument():
# Line 598  class Symbol(object): Line 685  class Symbol(object):
685            if isinstance(a,Symbol):            if isinstance(a,Symbol):
686               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
687            else:            else:
688                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
689                if len(s)>0:                if len(s)>0:
690                   out.append(numarray.zeros(s),numarray.Float64)                   out.append(numarray.zeros(s),numarray.Float64)
691                else:                else:
# Line 698  class Symbol(object): Line 785  class Symbol(object):
785         """         """
786         returns -self.         returns -self.
787    
788         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
789         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
790         """         """
791         return self*(-1.)         return self*(-1.)
# Line 707  class Symbol(object): Line 794  class Symbol(object):
794         """         """
795         returns +self.         returns +self.
796    
797         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
798         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
799         """         """
800         return self*(1.)         return self*(1.)
801    
802     def __abs__(self):     def __abs__(self):
803         """         """
804         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
805         """         """
806         return Abs_Symbol(self)         return Abs_Symbol(self)
807    
# Line 724  class Symbol(object): Line 811  class Symbol(object):
811    
812         @param other: object to be added to this object         @param other: object to be added to this object
813         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
814         @return:  a S{Symbol} representing the sum of this object and C{other}         @return:  a L{Symbol} representing the sum of this object and C{other}
815         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
816         """         """
817         return add(self,other)         return add(self,other)
# Line 735  class Symbol(object): Line 822  class Symbol(object):
822    
823         @param other: object this object is added to         @param other: object this object is added to
824         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
825         @return: a S{Symbol} representing the sum of C{other} and this object object         @return: a L{Symbol} representing the sum of C{other} and this object object
826         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
827         """         """
828         return add(other,self)         return add(other,self)
# Line 746  class Symbol(object): Line 833  class Symbol(object):
833    
834         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
835         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
836         @return: a S{Symbol} representing the difference of C{other} and this object         @return: a L{Symbol} representing the difference of C{other} and this object
837         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
838         """         """
839         return add(self,-other)         return add(self,-other)
# Line 757  class Symbol(object): Line 844  class Symbol(object):
844    
845         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
846         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
847         @return: a S{Symbol} representing the difference of this object and C{other}.         @return: a L{Symbol} representing the difference of this object and C{other}.
848         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
849         """         """
850         return add(-self,other)         return add(-self,other)
# Line 768  class Symbol(object): Line 855  class Symbol(object):
855    
856         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
857         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
858         @return: a S{Symbol} representing the product of the object and C{other}.         @return: a L{Symbol} representing the product of the object and C{other}.
859         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
860         """         """
861         return mult(self,other)         return mult(self,other)
# Line 779  class Symbol(object): Line 866  class Symbol(object):
866    
867         @param other: object this object is multiplied with         @param other: object this object is multiplied with
868         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
869         @return: a S{Symbol} representing the product of C{other} and the object.         @return: a L{Symbol} representing the product of C{other} and the object.
870         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
871         """         """
872         return mult(other,self)         return mult(other,self)
# Line 790  class Symbol(object): Line 877  class Symbol(object):
877    
878         @param other: object dividing this object         @param other: object dividing this object
879         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
880         @return: a S{Symbol} representing the quotient of this object and C{other}         @return: a L{Symbol} representing the quotient of this object and C{other}
881         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
882         """         """
883         return quotient(self,other)         return quotient(self,other)
# Line 801  class Symbol(object): Line 888  class Symbol(object):
888    
889         @param other: object dividing this object         @param other: object dividing this object
890         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
891         @return: a S{Symbol} representing the quotient of C{other} and this object         @return: a L{Symbol} representing the quotient of C{other} and this object
892         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
893         """         """
894         return quotient(other,self)         return quotient(other,self)
# Line 812  class Symbol(object): Line 899  class Symbol(object):
899    
900         @param other: exponent         @param other: exponent
901         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
902         @return: a S{Symbol} representing the power of this object to C{other}         @return: a L{Symbol} representing the power of this object to C{other}
903         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
904         """         """
905         return power(self,other)         return power(self,other)
# Line 823  class Symbol(object): Line 910  class Symbol(object):
910    
911         @param other: basis         @param other: basis
912         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
913         @return: a S{Symbol} representing the power of C{other} to this object         @return: a L{Symbol} representing the power of C{other} to this object
914         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
915         """         """
916         return power(other,self)         return power(other,self)
# Line 834  class Symbol(object): Line 921  class Symbol(object):
921    
922         @param index: defines a         @param index: defines a
923         @type index: C{slice} or C{int} or a C{tuple} of them         @type index: C{slice} or C{int} or a C{tuple} of them
924         @return: a S{Symbol} representing the slice defined by index         @return: a L{Symbol} representing the slice defined by index
925         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
926         """         """
927         return GetSlice_Symbol(self,index)         return GetSlice_Symbol(self,index)
# Line 844  class DependendSymbol(Symbol): Line 931  class DependendSymbol(Symbol):
931     DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.     DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.
932     Two DependendSymbol are equal if they have the same shape, the same arguments and one of them has an unspecified spatial dimension or the spatial dimension is identical       Two DependendSymbol are equal if they have the same shape, the same arguments and one of them has an unspecified spatial dimension or the spatial dimension is identical  
933        
934     Example:     Example::
935        
936     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
937     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
938     print u1==u2       print u1==u2
939     False       False
940        
941        but       but::
942    
943     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
944     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
945     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
946     print u1==u2, u1==u3       print u1==u2, u1==u3
947     True False       True False
948    
949     @note: DependendSymbol should be used as return value of functions with L{Symbol} arguments. This will allow the optimizer to remove redundant function calls.     @note: DependendSymbol should be used as return value of functions with L{Symbol} arguments. This will allow the optimizer to remove redundant function calls.
950     """     """
# Line 947  class GetSlice_Symbol(DependendSymbol): Line 1034  class GetSlice_Symbol(DependendSymbol):
1034        @type format: C{str}        @type format: C{str}
1035        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1036        @rtype: C{str}        @rtype: C{str}
1037        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1038        """        """
1039        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1040           return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])           return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
# Line 983  def log10(arg): Line 1070  def log10(arg):
1070    
1071     @param arg: argument     @param arg: argument
1072     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1073     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1074     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1075     """     """
1076     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1005  def wherePositive(arg): Line 1092  def wherePositive(arg):
1092    
1093     @param arg: argument     @param arg: argument
1094     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1095     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1096     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1097     """     """
1098     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1051  class WherePositive_Symbol(DependendSymb Line 1138  class WherePositive_Symbol(DependendSymb
1138        @type format: C{str}        @type format: C{str}
1139        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1140        @rtype: C{str}        @rtype: C{str}
1141        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1142        """        """
1143        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1144            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1087  def whereNegative(arg): Line 1174  def whereNegative(arg):
1174    
1175     @param arg: argument     @param arg: argument
1176     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1177     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1178     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1179     """     """
1180     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1133  class WhereNegative_Symbol(DependendSymb Line 1220  class WhereNegative_Symbol(DependendSymb
1220        @type format: C{str}        @type format: C{str}
1221        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1222        @rtype: C{str}        @rtype: C{str}
1223        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1224        """        """
1225        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1226            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1169  def whereNonNegative(arg): Line 1256  def whereNonNegative(arg):
1256    
1257     @param arg: argument     @param arg: argument
1258     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1259     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1260     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1261     """     """
1262     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1199  def whereNonPositive(arg): Line 1286  def whereNonPositive(arg):
1286    
1287     @param arg: argument     @param arg: argument
1288     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1289     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1290     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1291     """     """
1292     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1231  def whereZero(arg,tol=0.): Line 1318  def whereZero(arg,tol=0.):
1318     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1319     @param tol: tolerance. values with absolute value less then tol are accepted as zero.     @param tol: tolerance. values with absolute value less then tol are accepted as zero.
1320     @type tol: C{float}     @type tol: C{float}
1321     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1322     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1323     """     """
1324     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1239  def whereZero(arg,tol=0.): Line 1326  def whereZero(arg,tol=0.):
1326        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1327        return out        return out
1328     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1329        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1330     elif isinstance(arg,float):     elif isinstance(arg,float):
1331        if abs(arg)<=tol:        if abs(arg)<=tol:
1332          return 1.          return 1.
# Line 1280  class WhereZero_Symbol(DependendSymbol): Line 1364  class WhereZero_Symbol(DependendSymbol):
1364        @type format: C{str}        @type format: C{str}
1365        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1366        @rtype: C{str}        @rtype: C{str}
1367        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1368        """        """
1369        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1370           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1314  def whereNonZero(arg,tol=0.): Line 1398  def whereNonZero(arg,tol=0.):
1398    
1399     @param arg: argument     @param arg: argument
1400     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1401     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1402     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1403     """     """
1404     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1322  def whereNonZero(arg,tol=0.): Line 1406  def whereNonZero(arg,tol=0.):
1406        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1407        return out        return out
1408     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1409        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1410     elif isinstance(arg,float):     elif isinstance(arg,float):
1411        if abs(arg)>tol:        if abs(arg)>tol:
1412          return 1.          return 1.
# Line 1341  def whereNonZero(arg,tol=0.): Line 1422  def whereNonZero(arg,tol=0.):
1422     else:     else:
1423        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1424    
1425    def erf(arg):
1426       """
1427       returns erf of argument arg
1428    
1429       @param arg: argument
1430       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1431       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1432       @raises TypeError: if the type of the argument is not expected.
1433       """
1434       if isinstance(arg,escript.Data):
1435          return arg._erf()
1436       else:
1437          raise TypeError,"erf: Unknown argument type."
1438    
1439  def sin(arg):  def sin(arg):
1440     """     """
1441     returns sine of argument arg     returns sine of argument arg
1442    
1443     @param arg: argument     @param arg: argument
1444     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1445     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1446     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1447     """     """
1448     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1385  class Sin_Symbol(DependendSymbol): Line 1480  class Sin_Symbol(DependendSymbol):
1480        @type format: C{str}        @type format: C{str}
1481        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1482        @rtype: C{str}        @rtype: C{str}
1483        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1484        """        """
1485        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1486            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1437  def cos(arg): Line 1532  def cos(arg):
1532    
1533     @param arg: argument     @param arg: argument
1534     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1535     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1536     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1537     """     """
1538     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1475  class Cos_Symbol(DependendSymbol): Line 1570  class Cos_Symbol(DependendSymbol):
1570        @type format: C{str}        @type format: C{str}
1571        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1572        @rtype: C{str}        @rtype: C{str}
1573        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1574        """        """
1575        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1576            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1527  def tan(arg): Line 1622  def tan(arg):
1622    
1623     @param arg: argument     @param arg: argument
1624     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1625     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1626     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1627     """     """
1628     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1565  class Tan_Symbol(DependendSymbol): Line 1660  class Tan_Symbol(DependendSymbol):
1660        @type format: C{str}        @type format: C{str}
1661        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1662        @rtype: C{str}        @rtype: C{str}
1663        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1664        """        """
1665        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1666            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1617  def asin(arg): Line 1712  def asin(arg):
1712    
1713     @param arg: argument     @param arg: argument
1714     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1715     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1716     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1717     """     """
1718     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1655  class Asin_Symbol(DependendSymbol): Line 1750  class Asin_Symbol(DependendSymbol):
1750        @type format: C{str}        @type format: C{str}
1751        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1752        @rtype: C{str}        @rtype: C{str}
1753        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1754        """        """
1755        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1756            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1707  def acos(arg): Line 1802  def acos(arg):
1802    
1803     @param arg: argument     @param arg: argument
1804     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1805     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1806     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1807     """     """
1808     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1745  class Acos_Symbol(DependendSymbol): Line 1840  class Acos_Symbol(DependendSymbol):
1840        @type format: C{str}        @type format: C{str}
1841        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1842        @rtype: C{str}        @rtype: C{str}
1843        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1844        """        """
1845        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1846            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1797  def atan(arg): Line 1892  def atan(arg):
1892    
1893     @param arg: argument     @param arg: argument
1894     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1895     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1896     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1897     """     """
1898     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1835  class Atan_Symbol(DependendSymbol): Line 1930  class Atan_Symbol(DependendSymbol):
1930        @type format: C{str}        @type format: C{str}
1931        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1932        @rtype: C{str}        @rtype: C{str}
1933        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1934        """        """
1935        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1936            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1887  def sinh(arg): Line 1982  def sinh(arg):
1982    
1983     @param arg: argument     @param arg: argument
1984     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1985     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1986     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1987     """     """
1988     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1925  class Sinh_Symbol(DependendSymbol): Line 2020  class Sinh_Symbol(DependendSymbol):
2020        @type format: C{str}        @type format: C{str}
2021        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2022        @rtype: C{str}        @rtype: C{str}
2023        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2024        """        """
2025        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2026            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1977  def cosh(arg): Line 2072  def cosh(arg):
2072    
2073     @param arg: argument     @param arg: argument
2074     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2075     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2076     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2077     """     """
2078     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2015  class Cosh_Symbol(DependendSymbol): Line 2110  class Cosh_Symbol(DependendSymbol):
2110        @type format: C{str}        @type format: C{str}
2111        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2112        @rtype: C{str}        @rtype: C{str}
2113        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2114        """        """
2115        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2116            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2067  def tanh(arg): Line 2162  def tanh(arg):
2162    
2163     @param arg: argument     @param arg: argument
2164     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2165     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2166     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2167     """     """
2168     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2105  class Tanh_Symbol(DependendSymbol): Line 2200  class Tanh_Symbol(DependendSymbol):
2200        @type format: C{str}        @type format: C{str}
2201        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2202        @rtype: C{str}        @rtype: C{str}
2203        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2204        """        """
2205        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2206            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2157  def asinh(arg): Line 2252  def asinh(arg):
2252    
2253     @param arg: argument     @param arg: argument
2254     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2255     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2256     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2257     """     """
2258     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2195  class Asinh_Symbol(DependendSymbol): Line 2290  class Asinh_Symbol(DependendSymbol):
2290        @type format: C{str}        @type format: C{str}
2291        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2292        @rtype: C{str}        @rtype: C{str}
2293        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2294        """        """
2295        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2296            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2247  def acosh(arg): Line 2342  def acosh(arg):
2342    
2343     @param arg: argument     @param arg: argument
2344     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2345     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2346     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2347     """     """
2348     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2285  class Acosh_Symbol(DependendSymbol): Line 2380  class Acosh_Symbol(DependendSymbol):
2380        @type format: C{str}        @type format: C{str}
2381        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2382        @rtype: C{str}        @rtype: C{str}
2383        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2384        """        """
2385        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2386            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2337  def atanh(arg): Line 2432  def atanh(arg):
2432    
2433     @param arg: argument     @param arg: argument
2434     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2435     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2436     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2437     """     """
2438     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2375  class Atanh_Symbol(DependendSymbol): Line 2470  class Atanh_Symbol(DependendSymbol):
2470        @type format: C{str}        @type format: C{str}
2471        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2472        @rtype: C{str}        @rtype: C{str}
2473        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2474        """        """
2475        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2476            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2427  def exp(arg): Line 2522  def exp(arg):
2522    
2523     @param arg: argument     @param arg: argument
2524     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2525     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2526     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2527     """     """
2528     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2465  class Exp_Symbol(DependendSymbol): Line 2560  class Exp_Symbol(DependendSymbol):
2560        @type format: C{str}        @type format: C{str}
2561        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2562        @rtype: C{str}        @rtype: C{str}
2563        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2564        """        """
2565        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2566            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2517  def sqrt(arg): Line 2612  def sqrt(arg):
2612    
2613     @param arg: argument     @param arg: argument
2614     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2615     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2616     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2617     """     """
2618     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2555  class Sqrt_Symbol(DependendSymbol): Line 2650  class Sqrt_Symbol(DependendSymbol):
2650        @type format: C{str}        @type format: C{str}
2651        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2652        @rtype: C{str}        @rtype: C{str}
2653        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2654        """        """
2655        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2656            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2607  def log(arg): Line 2702  def log(arg):
2702    
2703     @param arg: argument     @param arg: argument
2704     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2705     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2706     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2707     """     """
2708     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2645  class Log_Symbol(DependendSymbol): Line 2740  class Log_Symbol(DependendSymbol):
2740        @type format: C{str}        @type format: C{str}
2741        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2742        @rtype: C{str}        @rtype: C{str}
2743        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2744        """        """
2745        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2746            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2697  def sign(arg): Line 2792  def sign(arg):
2792    
2793     @param arg: argument     @param arg: argument
2794     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2795     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2796     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2797     """     """
2798     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2745  class Abs_Symbol(DependendSymbol): Line 2840  class Abs_Symbol(DependendSymbol):
2840        @type format: C{str}        @type format: C{str}
2841        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2842        @rtype: C{str}        @rtype: C{str}
2843        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2844        """        """
2845        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2846            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2797  def minval(arg): Line 2892  def minval(arg):
2892    
2893     @param arg: argument     @param arg: argument
2894     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2895     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2896     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2897     """     """
2898     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2838  class Minval_Symbol(DependendSymbol): Line 2933  class Minval_Symbol(DependendSymbol):
2933        @type format: C{str}        @type format: C{str}
2934        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2935        @rtype: C{str}        @rtype: C{str}
2936        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2937        """        """
2938        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2939            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2874  def maxval(arg): Line 2969  def maxval(arg):
2969    
2970     @param arg: argument     @param arg: argument
2971     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2972     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2973     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2974     """     """
2975     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2915  class Maxval_Symbol(DependendSymbol): Line 3010  class Maxval_Symbol(DependendSymbol):
3010        @type format: C{str}        @type format: C{str}
3011        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3012        @rtype: C{str}        @rtype: C{str}
3013        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3014        """        """
3015        if isinstance(argstrs,list):        if isinstance(argstrs,list):
3016            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2951  def length(arg): Line 3046  def length(arg):
3046    
3047     @param arg: argument     @param arg: argument
3048     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3049     @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
3050     """     """
3051     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
3052    
# Line 2961  def trace(arg,axis_offset=0): Line 3056  def trace(arg,axis_offset=0):
3056    
3057     @param arg: argument     @param arg: argument
3058     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3059     @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
3060                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3061     @type axis_offset: C{int}     @type axis_offset: C{int}
3062     @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.
3063     @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 2970  def trace(arg,axis_offset=0): Line 3065  def trace(arg,axis_offset=0):
3065     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
3066        sh=arg.shape        sh=arg.shape
3067        if len(sh)<2:        if len(sh)<2:
3068          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3069        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
3070          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
3071        s1=1        s1=1
3072        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
3073        s2=1        s2=1
3074        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3075        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
3076          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)
3077        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))
3078        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
3079        for i1 in range(s1):        for i1 in range(s1):
# Line 2987  def trace(arg,axis_offset=0): Line 3082  def trace(arg,axis_offset=0):
3082        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3083        return out        return out
3084     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3085        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
3086            raise ValueError,"rank of argument must be greater than 1"
3087          if axis_offset<0 or axis_offset>arg.getRank()-2:
3088            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3089          s=list(arg.getShape())        
3090          if not s[axis_offset] == s[axis_offset+1]:
3091            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3092          return arg._trace(axis_offset)
3093     elif isinstance(arg,float):     elif isinstance(arg,float):
3094        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
3095     elif isinstance(arg,int):     elif isinstance(arg,int):
3096        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
3097     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3098        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
3099     else:     else:
3100        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
3101    
 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  
3102  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
3103     """     """
3104     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 3045  class Trace_Symbol(DependendSymbol): Line 3108  class Trace_Symbol(DependendSymbol):
3108        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3109        @param arg: argument of function        @param arg: argument of function
3110        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3111        @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
3112                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3113        @type axis_offset: C{int}        @type axis_offset: C{int}
3114        """        """
3115        if arg.getRank()<2:        if arg.getRank()<2:
3116          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3117        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3118          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
3119        s=list(arg.getShape())                s=list(arg.getShape())        
3120        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3121          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)
3122        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())
3123    
3124     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 3068  class Trace_Symbol(DependendSymbol): Line 3131  class Trace_Symbol(DependendSymbol):
3131        @type format: C{str}        @type format: C{str}
3132        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3133        @rtype: C{str}        @rtype: C{str}
3134        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3135        """        """
3136        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3137           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3112  class Trace_Symbol(DependendSymbol): Line 3175  class Trace_Symbol(DependendSymbol):
3175    
3176  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3177     """     """
3178     returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.     returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3179    
3180     @param arg: argument     @param arg: argument
3181     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3182     @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.     @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3183                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3184     @type axis_offset: C{int}     @type axis_offset: C{int}
3185     @return: transpose of arg     @return: transpose of arg
3186     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
# Line 3126  def transpose(arg,axis_offset=None): Line 3189  def transpose(arg,axis_offset=None):
3189        if axis_offset==None: axis_offset=int(arg.rank/2)        if axis_offset==None: axis_offset=int(arg.rank/2)
3190        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3191     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3192        if axis_offset==None: axis_offset=int(arg.getRank()/2)        r=arg.getRank()
3193        return escript_transpose(arg,axis_offset)        if axis_offset==None: axis_offset=int(r/2)
3194          if axis_offset<0 or axis_offset>r:
3195            raise ValueError,"axis_offset must be between 0 and %s"%r
3196          return arg._transpose(axis_offset)
3197     elif isinstance(arg,float):     elif isinstance(arg,float):
3198        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3199          raise ValueError,"transpose: axis_offset must be 0 for float argument"          raise ValueError,"axis_offset must be 0 for float argument"
3200        return arg        return arg
3201     elif isinstance(arg,int):     elif isinstance(arg,int):
3202        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3203          raise ValueError,"transpose: axis_offset must be 0 for int argument"          raise ValueError,"axis_offset must be 0 for int argument"
3204        return float(arg)        return float(arg)
3205     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3206        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3207        return Transpose_Symbol(arg,axis_offset)        return Transpose_Symbol(arg,axis_offset)
3208     else:     else:
3209        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"Unknown argument type."
3210    
 def escript_transpose(arg,axis_offset): # this should be escript._transpose  
       "arg si a Data objects!!!"  
       r=arg.getRank()  
       if axis_offset<0 or axis_offset>r:  
         raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r  
       s=arg.getShape()  
       s_out=s[axis_offset:]+s[:axis_offset]  
       out=escript.Data(0.,s_out,arg.getFunctionSpace())  
       if r==4:  
          if axis_offset==1:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i3,i0,i1,i2]  
          elif axis_offset==2:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i2,i3,i0,i1]  
          elif axis_offset==3:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i1,i2,i3,i0]  
          else:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i0,i1,i2,i3]  
       elif r==3:  
          if axis_offset==1:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                          out[i0,i1,i2]=arg[i2,i0,i1]  
          elif axis_offset==2:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                          out[i0,i1,i2]=arg[i1,i2,i0]  
          else:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                          out[i0,i1,i2]=arg[i0,i1,i2]  
       elif r==2:  
          if axis_offset==1:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                          out[i0,i1]=arg[i1,i0]  
          else:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                          out[i0,i1]=arg[i0,i1]  
       elif r==1:  
           for i0 in range(s_out[0]):  
                out[i0]=arg[i0]  
       elif r==0:  
              out=arg+0.  
       return out  
3211  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3212     """     """
3213     L{Symbol} representing the result of the transpose function     L{Symbol} representing the result of the transpose function
# Line 3216  class Transpose_Symbol(DependendSymbol): Line 3218  class Transpose_Symbol(DependendSymbol):
3218    
3219        @param arg: argument of function        @param arg: argument of function
3220        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3221         @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.        @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3222                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3223        @type axis_offset: C{int}        @type axis_offset: C{int}
3224        """        """
3225        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3226        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3227          raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3228        s=arg.getShape()        s=arg.getShape()
3229        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3230    
# Line 3236  class Transpose_Symbol(DependendSymbol): Line 3238  class Transpose_Symbol(DependendSymbol):
3238        @type format: C{str}        @type format: C{str}
3239        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3240        @rtype: C{str}        @rtype: C{str}
3241        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3242        """        """
3243        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3244           return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3277  class Transpose_Symbol(DependendSymbol): Line 3279  class Transpose_Symbol(DependendSymbol):
3279           return identity(self.getShape())           return identity(self.getShape())
3280        else:        else:
3281           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3282    
3283    def swap_axes(arg,axis0=0,axis1=1):
3284       """
3285       returns the swap of arg by swaping the components axis0 and axis1
3286    
3287       @param arg: argument
3288       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3289       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3290       @type axis0: C{int}
3291       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3292       @type axis1: C{int}
3293       @return: C{arg} with swaped components
3294       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3295       """
3296       if axis0 > axis1:
3297          axis0,axis1=axis1,axis0
3298       if isinstance(arg,numarray.NumArray):
3299          return numarray.swapaxes(arg,axis0,axis1)
3300       elif isinstance(arg,escript.Data):
3301          return arg._swap_axes(axis0,axis1)
3302       elif isinstance(arg,float):
3303          raise TyepError,"float argument is not supported."
3304       elif isinstance(arg,int):
3305          raise TyepError,"int argument is not supported."
3306       elif isinstance(arg,Symbol):
3307          return SwapAxes_Symbol(arg,axis0,axis1)
3308       else:
3309          raise TypeError,"Unknown argument type."
3310    
3311    class SwapAxes_Symbol(DependendSymbol):
3312       """
3313       L{Symbol} representing the result of the swap function
3314       """
3315       def __init__(self,arg,axis0=0,axis1=1):
3316          """
3317          initialization of swap L{Symbol} with argument arg
3318    
3319          @param arg: argument
3320          @type arg: L{Symbol}.
3321          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3322          @type axis0: C{int}
3323          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3324          @type axis1: C{int}
3325          """
3326          if arg.getRank()<2:
3327             raise ValueError,"argument must have at least rank 2."
3328          if axis0<0 or axis0>arg.getRank()-1:
3329             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3330          if axis1<0 or axis1>arg.getRank()-1:
3331             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3332          if axis0 == axis1:
3333             raise ValueError,"axis indices must be different."
3334          if axis0 > axis1:
3335             axis0,axis1=axis1,axis0
3336          s=arg.getShape()
3337          s_out=[]
3338          for i in range(len(s)):
3339             if i == axis0:
3340                s_out.append(s[axis1])
3341             elif i == axis1:
3342                s_out.append(s[axis0])
3343             else:
3344                s_out.append(s[i])
3345          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3346    
3347       def getMyCode(self,argstrs,format="escript"):
3348          """
3349          returns a program code that can be used to evaluate the symbol.
3350    
3351          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3352          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3353          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3354          @type format: C{str}
3355          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3356          @rtype: C{str}
3357          @raise NotImplementedError: if the requested format is not available
3358          """
3359          if format=="escript" or format=="str"  or format=="text":
3360             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3361          else:
3362             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3363    
3364       def substitute(self,argvals):
3365          """
3366          assigns new values to symbols in the definition of the symbol.
3367          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3368    
3369          @param argvals: new values assigned to symbols
3370          @type argvals: C{dict} with keywords of type L{Symbol}.
3371          @return: result of the substitution process. Operations are executed as much as possible.
3372          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3373          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3374          """
3375          if argvals.has_key(self):
3376             arg=argvals[self]
3377             if self.isAppropriateValue(arg):
3378                return arg
3379             else:
3380                raise TypeError,"%s: new value is not appropriate."%str(self)
3381          else:
3382             arg=self.getSubstitutedArguments(argvals)
3383             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3384    
3385       def diff(self,arg):
3386          """
3387          differential of this object
3388    
3389          @param arg: the derivative is calculated with respect to arg
3390          @type arg: L{escript.Symbol}
3391          @return: derivative with respect to C{arg}
3392          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3393          """
3394          if arg==self:
3395             return identity(self.getShape())
3396          else:
3397             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3398    
3399  def symmetric(arg):  def symmetric(arg):
3400      """      """
3401      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
# Line 3289  def symmetric(arg): Line 3408  def symmetric(arg):
3408      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3409        if arg.rank==2:        if arg.rank==2:
3410          if not (arg.shape[0]==arg.shape[1]):          if not (arg.shape[0]==arg.shape[1]):
3411             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3412        elif arg.rank==4:        elif arg.rank==4:
3413          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3414             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3415        else:        else:
3416          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3417        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3418      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3419        return escript_symmetric(arg)        if arg.getRank()==2:
3420            if not (arg.getShape()[0]==arg.getShape()[1]):
3421               raise ValueError,"argument must be square."
3422            return arg._symmetric()
3423          elif arg.getRank()==4:
3424            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3425               raise ValueError,"argument must be square."
3426            return arg._symmetric()
3427          else:
3428            raise ValueError,"rank 2 or 4 is required."
3429      elif isinstance(arg,float):      elif isinstance(arg,float):
3430        return arg        return arg
3431      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3305  def symmetric(arg): Line 3433  def symmetric(arg):
3433      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3434        if arg.getRank()==2:        if arg.getRank()==2:
3435          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3436             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3437        elif arg.getRank()==4:        elif arg.getRank()==4:
3438          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3439             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3440        else:        else:
3441          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3442        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3443      else:      else:
3444        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3445    
 def escript_symmetric(arg): # this should be implemented in c++  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               out[i0,i1]=(arg[i0,i1]+arg[i1,i0])/2.  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               for i2 in range(arg.getShape()[2]):  
                  for i3 in range(arg.getShape()[3]):  
                      out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]+arg[i2,i3,i0,i1])/2.  
       else:  
         raise ValueError,"escript_symmetric: rank 2 or 4 is required."  
       return out  
   
3446  def nonsymmetric(arg):  def nonsymmetric(arg):
3447      """      """
3448      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
# Line 3356  def nonsymmetric(arg): Line 3463  def nonsymmetric(arg):
3463          raise ValueError,"nonsymmetric: rank 2 or 4 is required."          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3464        return (arg-transpose(arg))/2        return (arg-transpose(arg))/2
3465      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3466        return escript_nonsymmetric(arg)        if arg.getRank()==2:
3467            if not (arg.getShape()[0]==arg.getShape()[1]):
3468               raise ValueError,"argument must be square."
3469            return arg._nonsymmetric()
3470          elif arg.getRank()==4:
3471            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3472               raise ValueError,"argument must be square."
3473            return arg._nonsymmetric()
3474          else:
3475            raise ValueError,"rank 2 or 4 is required."
3476      elif isinstance(arg,float):      elif isinstance(arg,float):
3477        return arg        return arg
3478      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3374  def nonsymmetric(arg): Line 3490  def nonsymmetric(arg):
3490      else:      else:
3491        raise TypeError,"nonsymmetric: Unknown argument type."        raise TypeError,"nonsymmetric: Unknown argument type."
3492    
 def escript_nonsymmetric(arg): # this should be implemented in c++  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               out[i0,i1]=(arg[i0,i1]-arg[i1,i0])/2.  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               for i2 in range(arg.getShape()[2]):  
                  for i3 in range(arg.getShape()[3]):  
                      out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]-arg[i2,i3,i0,i1])/2.  
       else:  
         raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."  
       return out  
   
   
3493  def inverse(arg):  def inverse(arg):
3494      """      """
3495      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3496    
3497      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3498      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3499      @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])      @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3500      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3501      @remark: for L{escript.Data} objects the dimension is restricted to 3.      @note: for L{escript.Data} objects the dimension is restricted to 3.
3502      """      """
3503        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3504      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3505        return numarray.linear_algebra.inverse(arg)        return numarray.linear_algebra.inverse(arg)
3506      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
# Line 3498  class Inverse_Symbol(DependendSymbol): Line 3593  class Inverse_Symbol(DependendSymbol):
3593        @type format: C{str}        @type format: C{str}
3594        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3595        @rtype: C{str}        @rtype: C{str}
3596        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3597        """        """
3598        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3599           return "inverse(%s)"%argstrs[0]           return "inverse(%s)"%argstrs[0]
# Line 3538  class Inverse_Symbol(DependendSymbol): Line 3633  class Inverse_Symbol(DependendSymbol):
3633        if arg==self:        if arg==self:
3634           return identity(self.getShape())           return identity(self.getShape())
3635        else:        else:
3636           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3637    
3638  def eigenvalues(arg):  def eigenvalues(arg):
3639      """      """
# Line 3549  def eigenvalues(arg): Line 3644  def eigenvalues(arg):
3644      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3645      @return: the eigenvalues in increasing order.      @return: the eigenvalues in increasing order.
3646      @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.
3647      @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.      @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3648      """      """
3649      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3650        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
# Line 3610  def eigenvalues(arg): Line 3705  def eigenvalues(arg):
3705    
3706  def eigenvalues_and_eigenvectors(arg):  def eigenvalues_and_eigenvectors(arg):
3707      """      """
3708      returns the eigenvalues of the square matrix arg.      returns the eigenvalues and eigenvectors of the square matrix arg.
3709    
3710      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3711                  arg must be symmetric, ie. transpose(arg)==arg (this is not checked).                  arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3712      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{escript.Data}
3713      @return: the eigenvalues in increasing order.      @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3714      @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.               eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3715      @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.               the eigenvector coresponding to the i-th eigenvalue.
3716        @rtype: L{tuple} of L{escript.Data}.
3717        @note: The dimension is restricted to 3.
3718      """      """
3719      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3720        raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"        raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
# Line 3674  class Add_Symbol(DependendSymbol): Line 3771  class Add_Symbol(DependendSymbol):
3771         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3772         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3773         """         """
3774         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3775         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3776         if not sh0==sh1:         if not sh0==sh1:
3777            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3778         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3690  class Add_Symbol(DependendSymbol): Line 3787  class Add_Symbol(DependendSymbol):
3787        @type format: C{str}        @type format: C{str}
3788        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.        @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3789        @rtype: C{str}        @rtype: C{str}
3790        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3791        """        """
3792        if format=="str" or format=="text":        if format=="str" or format=="text":
3793           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 3749  def mult(arg0,arg1): Line 3846  def mult(arg0,arg1):
3846         """         """
3847         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3848         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3849            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3850         else:         else:
3851            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3852                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3773  class Mult_Symbol(DependendSymbol): Line 3870  class Mult_Symbol(DependendSymbol):
3870         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3871         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3872         """         """
3873         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3874         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3875         if not sh0==sh1:         if not sh0==sh1:
3876            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3877         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 3789  class Mult_Symbol(DependendSymbol): Line 3886  class Mult_Symbol(DependendSymbol):
3886        @type format: C{str}        @type format: C{str}
3887        @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.
3888        @rtype: C{str}        @rtype: C{str}
3889        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3890        """        """
3891        if format=="str" or format=="text":        if format=="str" or format=="text":
3892           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3849  def quotient(arg0,arg1): Line 3946  def quotient(arg0,arg1):
3946         """         """
3947         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3948         if testForZero(args[0]):         if testForZero(args[0]):
3949            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3950         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3951            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3952               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3878  class Quotient_Symbol(DependendSymbol): Line 3975  class Quotient_Symbol(DependendSymbol):
3975         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3976         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3977         """         """
3978         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3979         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3980         if not sh0==sh1:         if not sh0==sh1:
3981            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3982         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 3894  class Quotient_Symbol(DependendSymbol): Line 3991  class Quotient_Symbol(DependendSymbol):
3991        @type format: C{str}        @type format: C{str}
3992        @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.
3993        @rtype: C{str}        @rtype: C{str}
3994        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3995        """        """
3996        if format=="str" or format=="text":        if format=="str" or format=="text":
3997           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3955  def power(arg0,arg1): Line 4052  def power(arg0,arg1):
4052         """         """
4053         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4054         if testForZero(args[0]):         if testForZero(args[0]):
4055            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4056         elif testForZero(args[1]):         elif testForZero(args[1]):
4057            return numarray.ones(pokeShape(args[1]),numarray.Float64)            return numarray.ones(getShape(args[1]),numarray.Float64)
4058         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4059            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4060         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 3980  class Power_Symbol(DependendSymbol): Line 4077  class Power_Symbol(DependendSymbol):
4077         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4078         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4079         """         """
4080         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4081         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4082         if not sh0==sh1:         if not sh0==sh1:
4083            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4084         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3998  class Power_Symbol(DependendSymbol): Line 4095  class Power_Symbol(DependendSymbol):
4095        @type format: C{str}        @type format: C{str}
4096        @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.
4097        @rtype: C{str}        @rtype: C{str}
4098        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4099        """        """
4100        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4101           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 4080  def minimum(*args): Line 4177  def minimum(*args):
4177            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4178      return out      return out
4179    
4180  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4181      """      """
4182      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4183    
4184      @param arg: argument      @param arg: argument
4185      @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}
4186      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4187      @type arg: C{float}      @type minval: C{float} or C{None}
4188      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4189      @type arg: C{float}      @type maxval: C{float} or C{None}
4190      @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
4191               less then maxval are unchanged.               less then maxval are unchanged.
4192      @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
4193      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4194      """      """
4195      if minval>maxval:      if not minval==None and not maxval==None:
4196         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4197      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4198        if minval == None:
4199            tmp=arg
4200        else:
4201            tmp=maximum(minval,arg)
4202        if maxval == None:
4203            return tmp
4204        else:
4205            return minimum(tmp,maxval)
4206    
4207        
4208  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 4114  def inner(arg0,arg1): Line 4219  def inner(arg0,arg1):
4219      @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}
4220      @param arg1: second argument      @param arg1: second argument
4221      @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}
4222      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4223      @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
4224      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4225      """      """
4226      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4227      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4228      if not sh0==sh1:      if not sh0==sh1:
4229          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4230      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4231    
4232    def outer(arg0,arg1):
4233        """
4234        the outer product of the two argument:
4235    
4236        out[t,s]=arg0[t]*arg1[s]
4237    
4238        where
4239    
4240            - s runs through arg0.Shape
4241            - t runs through arg1.Shape
4242    
4243        @param arg0: first argument
4244        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4245        @param arg1: second argument
4246        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4247        @return: the outer product of arg0 and arg1 at each data point
4248        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4249        """
4250        return generalTensorProduct(arg0,arg1,axis_offset=0)
4251    
4252  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4253      """      """
4254        see L{matrix_mult}
4255        """
4256        return matrix_mult(arg0,arg1)
4257    
4258    def matrix_mult(arg0,arg1):
4259        """
4260      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4261    
4262      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4263    
4264            or      or
4265    
4266      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4267    
4268      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.
4269    
4270      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4271      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4144  def matrixmult(arg0,arg1): Line 4275  def matrixmult(arg0,arg1):
4275      @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
4276      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4277      """      """
4278      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4279      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4280      if not len(sh0)==2 :      if not len(sh0)==2 :
4281          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4282      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4283          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4284      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4285    
4286  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4287      """      """
4288      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  
4289      """      """
4290      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
4291    
4292    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4293      """      """
4294      the tensor product of the two argument:      the tensor product of the two argument:
   
4295            
4296      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4297    
4298      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4299    
4300                   or      or
4301    
4302      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4303    
# Line 4189  def tensormult(arg0,arg1): Line 4306  def tensormult(arg0,arg1):
4306    
4307      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]
4308                                
4309                   or      or
4310    
4311      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]
4312    
4313                   or      or
4314    
4315      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]
4316    
4317      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  
4318      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.
4319    
4320      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4321      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4207  def tensormult(arg0,arg1): Line 4324  def tensormult(arg0,arg1):
4324      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4325      @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
4326      """      """
4327      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4328      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4329      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4330         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4331      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):
4332         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4333      else:      else:
4334          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4335    
4336  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4337      """      """
# Line 4222  def generalTensorProduct(arg0,arg1,axis_ Line 4339  def generalTensorProduct(arg0,arg1,axis_
4339    
4340      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4341    
4342      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:]  
4343    
4344      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]
4345      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4346            - t runs through arg1.Shape[axis_offset:]
4347    
4348      @param arg0: first argument      @param arg0: first argument
4349      @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}
4350      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4351      @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}
4352      @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.
4353      @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 4244  def generalTensorProduct(arg0,arg1,axis_ Line 4360  def generalTensorProduct(arg0,arg1,axis_
4360             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4361         else:         else:
4362             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4363                 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)
4364             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4365             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4366             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 4270  def generalTensorProduct(arg0,arg1,axis_ Line 4386  def generalTensorProduct(arg0,arg1,axis_
4386                                    
4387  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4388     """     """
4389     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4390     """     """
4391     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4392         """         """
4393         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4394    
4395         @param arg0: numerator         @param arg0: first argument
4396         @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}.
4397         @param arg1: denominator         @param arg1: second argument
4398         @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}.
4399         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4400         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4401         """         """
4402         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4403         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4404         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4405         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4406         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 4303  class GeneralTensorProduct_Symbol(Depend Line 4419  class GeneralTensorProduct_Symbol(Depend
4419        @type format: C{str}        @type format: C{str}
4420        @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.
4421        @rtype: C{str}        @rtype: C{str}
4422        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4423        """        """
4424        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4425           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 4331  class GeneralTensorProduct_Symbol(Depend Line 4447  class GeneralTensorProduct_Symbol(Depend
4447           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4448           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4449    
4450  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4451      "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!!!"
4452      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4453      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
4454      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  def transposed_matrix_mult(arg0,arg1):
4455      shape10=arg1.getShape()[:axis_offset]      """
4456      shape1=arg1.getShape()[axis_offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4457      if not shape01==shape10:  
4458          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]
   
     # whatr function space should be used? (this here is not good!)  
     fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
     # create return value:  
     out=escript.Data(0.,tuple(shape0+shape1),fs)  
     #  
     s0=[[]]  
     for k in shape0:  
           s=[]  
           for j in s0:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s0=s  
     s1=[[]]  
     for k in shape1:  
           s=[]  
           for j in s1:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s1=s  
     s01=[[]]  
     for k in shape01:  
           s=[]  
           for j in s01:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s01=s  
   
     for i0 in s0:  
        for i1 in s1:  
          s=escript.Scalar(0.,fs)  
          for i01 in s01:  
             s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))  
          out.__setitem__(tuple(i0+i1),s)  
     return out  
4459    
4460        or
4461    
4462        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4463    
4464        The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4465    
4466        The first dimension of arg0 and arg1 must match.
4467    
4468        @param arg0: first argument of rank 2
4469        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4470        @param arg1: second argument of at least rank 1
4471        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4472        @return: the product of the transposed of arg0 and arg1 at each data point
4473        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4474        @raise ValueError: if the shapes of the arguments are not appropriate
4475        """
4476        sh0=getShape(arg0)
4477        sh1=getShape(arg1)
4478        if not len(sh0)==2 :
4479            raise ValueError,"first argument must have rank 2"
4480        if not len(sh1)==2 and not len(sh1)==1:
4481            raise ValueError,"second argument must have rank 1 or 2"
4482        return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4483    
4484    def transposed_tensor_mult(arg0,arg1):
4485        """
4486        the tensor product of the transposed of the first and the second argument
4487        
4488        for arg0 of rank 2 this is
4489    
4490        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4491    
4492        or
4493    
4494        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4495    
4496      
4497        and for arg0 of rank 4 this is
4498    
4499        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4500                  
4501        or
4502    
4503        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4504    
4505        or
4506    
4507        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4508    
4509        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4510        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4511    
4512        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4513    
4514        @param arg0: first argument of rank 2 or 4
4515        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4516        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4517        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4518        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4519        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4520        """
4521        sh0=getShape(arg0)
4522        sh1=getShape(arg1)
4523        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4524           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4525        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4526           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4527        else:
4528            raise ValueError,"first argument must have rank 2 or 4"
4529    
4530    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4531        """
4532        generalized tensor product of transposed of arg0 and arg1:
4533    
4534        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4535    
4536        where
4537    
4538            - s runs through arg0.Shape[axis_offset:]
4539            - r runs trough arg0.Shape[:axis_offset]
4540            - t runs through arg1.Shape[axis_offset:]
4541    
4542        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4543        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4544    
4545        @param arg0: first argument
4546        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4547        @param arg1: second argument
4548        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4549        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4550        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4551        """
4552        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4553        arg0,arg1=matchType(arg0,arg1)
4554        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4555        if isinstance(arg0,numarray.NumArray):
4556           if isinstance(arg1,Symbol):
4557               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4558           else:
4559               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4560                   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)
4561               arg0_c=arg0.copy()
4562               arg1_c=arg1.copy()
4563               sh0,sh1=arg0.shape,arg1.shape
4564               d0,d1,d01=1,1,1
4565               for i in sh0[axis_offset:]: d0*=i
4566               for i in sh1[axis_offset:]: d1*=i
4567               for i in sh0[:axis_offset]: d01*=i
4568               arg0_c.resize((d01,d0))
4569               arg1_c.resize((d01,d1))
4570               out=numarray.zeros((d0,d1),numarray.Float64)
4571               for i0 in range(d0):
4572                        for i1 in range(d1):
4573                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4574               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4575               return out
4576        elif isinstance(arg0,escript.Data):
4577           if isinstance(arg1,Symbol):
4578               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4579           else:
4580               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4581        else:      
4582           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4583                    
4584    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4585       """
4586       Symbol representing the general tensor product of the transposed of arg0 and arg1
4587       """
4588       def __init__(self,arg0,arg1,axis_offset=0):
4589           """
4590           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4591    
4592           @param arg0: first argument
4593           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4594           @param arg1: second argument
4595           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4596           @raise ValueError: inconsistent dimensions of arguments.
4597           @note: if both arguments have a spatial dimension, they must equal.
4598           """
4599           sh_arg0=getShape(arg0)
4600           sh_arg1=getShape(arg1)
4601           sh01=sh_arg0[:axis_offset]
4602           sh10=sh_arg1[:axis_offset]
4603           sh0=sh_arg0[axis_offset:]
4604           sh1=sh_arg1[axis_offset:]
4605           if not sh01==sh10:
4606               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)
4607           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4608    
4609       def getMyCode(self,argstrs,format="escript"):
4610          """
4611          returns a program code that can be used to evaluate the symbol.
4612    
4613          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4614          @type argstrs: C{list} of length 2 of C{str}.
4615          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4616          @type format: C{str}
4617          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4618          @rtype: C{str}
4619          @raise NotImplementedError: if the requested format is not available
4620          """
4621          if format=="escript" or format=="str" or format=="text":
4622             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4623          else:
4624             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4625    
4626       def substitute(self,argvals):
4627          """
4628          assigns new values to symbols in the definition of the symbol.
4629          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4630    
4631          @param argvals: new values assigned to symbols
4632          @type argvals: C{dict} with keywords of type L{Symbol}.
4633          @return: result of the substitution process. Operations are executed as much as possible.
4634          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4635          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4636          """
4637          if argvals.has_key(self):
4638             arg=argvals[self]
4639             if self.isAppropriateValue(arg):
4640                return arg
4641             else:
4642                raise TypeError,"%s: new value is not appropriate."%str(self)
4643          else:
4644             args=self.getSubstitutedArguments(argvals)
4645             return generalTransposedTensorProduct(args[0],args[1],args[2])
4646    
4647    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4648        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4649        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4650    
4651    def matrix_transposed_mult(arg0,arg1):
4652        """
4653        matrix-transposed(matrix) product of the two argument:
4654    
4655        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4656    
4657        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4658    
4659        The last dimensions of arg0 and arg1 must match.
4660    
4661        @param arg0: first argument of rank 2
4662        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4663        @param arg1: second argument of rank 2
4664        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4665        @return: the product of arg0 and the transposed of arg1 at each data point
4666        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4667        @raise ValueError: if the shapes of the arguments are not appropriate
4668        """
4669        sh0=getShape(arg0)
4670        sh1=getShape(arg1)
4671        if not len(sh0)==2 :
4672            raise ValueError,"first argument must have rank 2"
4673        if not len(sh1)==2 and not len(sh1)==1:
4674            raise ValueError,"second argument must have rank 1 or 2"
4675        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4676    
4677    def tensor_transposed_mult(arg0,arg1):
4678        """
4679        the tensor product of the first and the transpose of the second argument
4680        
4681        for arg0 of rank 2 this is
4682    
4683        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4684    
4685        and for arg0 of rank 4 this is
4686    
4687        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4688                  
4689        or
4690    
4691        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4692    
4693        In the first case the the second dimension of arg0 and arg1 must match and  
4694        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4695    
4696        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4697    
4698        @param arg0: first argument of rank 2 or 4
4699        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4700        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4701        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4702        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4703        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4704        """
4705        sh0=getShape(arg0)
4706        sh1=getShape(arg1)
4707        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4708           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4709        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4710           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4711        else:
4712            raise ValueError,"first argument must have rank 2 or 4"
4713    
4714    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4715        """
4716        generalized tensor product of transposed of arg0 and arg1:
4717    
4718        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4719    
4720        where
4721    
4722            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4723            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4724            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4725    
4726        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4727        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4728    
4729        @param arg0: first argument
4730        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4731        @param arg1: second argument
4732        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4733        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4734        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4735        """
4736        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4737        arg0,arg1=matchType(arg0,arg1)
4738        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4739        if isinstance(arg0,numarray.NumArray):
4740           if isinstance(arg1,Symbol):
4741               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4742           else:
4743               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4744                   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)
4745               arg0_c=arg0.copy()
4746               arg1_c=arg1.copy()
4747               sh0,sh1=arg0.shape,arg1.shape
4748               d0,d1,d01=1,1,1
4749               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4750               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4751               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4752               arg0_c.resize((d0,d01))
4753               arg1_c.resize((d1,d01))
4754               out=numarray.zeros((d0,d1),numarray.Float64)
4755               for i0 in range(d0):
4756                        for i1 in range(d1):
4757                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4758               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4759               return out
4760        elif isinstance(arg0,escript.Data):
4761           if isinstance(arg1,Symbol):
4762               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4763           else:
4764               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4765        else:      
4766           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4767                    
4768    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4769       """
4770       Symbol representing the general tensor product of arg0 and the transpose of arg1
4771       """
4772       def __init__(self,arg0,arg1,axis_offset=0):
4773           """
4774           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4775    
4776           @param arg0: first argument
4777           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4778           @param arg1: second argument
4779           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4780           @raise ValueError: inconsistent dimensions of arguments.
4781           @note: if both arguments have a spatial dimension, they must equal.
4782           """
4783           sh_arg0=getShape(arg0)
4784           sh_arg1=getShape(arg1)
4785           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4786           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4787           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4788           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4789           if not sh01==sh10:
4790               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)
4791           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4792    
4793       def getMyCode(self,argstrs,format="escript"):
4794          """
4795          returns a program code that can be used to evaluate the symbol.
4796    
4797          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4798          @type argstrs: C{list} of length 2 of C{str}.
4799          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4800          @type format: C{str}
4801          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4802          @rtype: C{str}
4803          @raise NotImplementedError: if the requested format is not available
4804          """
4805          if format=="escript" or format=="str" or format=="text":
4806             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4807          else:
4808             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4809    
4810       def substitute(self,argvals):
4811          """
4812          assigns new values to symbols in the definition of the symbol.
4813          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4814    
4815          @param argvals: new values assigned to symbols
4816          @type argvals: C{dict} with keywords of type L{Symbol}.
4817          @return: result of the substitution process. Operations are executed as much as possible.
4818          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4819          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4820          """
4821          if argvals.has_key(self):
4822             arg=argvals[self]
4823             if self.isAppropriateValue(arg):
4824                return arg
4825             else:
4826                raise TypeError,"%s: new value is not appropriate."%str(self)
4827          else:
4828             args=self.getSubstitutedArguments(argvals)
4829             return generalTensorTransposedProduct(args[0],args[1],args[2])
4830    
4831    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4832        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4833        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4834    
4835  #=========================================================  #=========================================================
4836  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
# Line 4394  def grad(arg,where=None): Line 4852  def grad(arg,where=None):
4852                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4853      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4854      @return: gradient of arg.      @return: gradient of arg.
4855      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4856      """      """
4857      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4858         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 4434  class Grad_Symbol(DependendSymbol): Line 4892  class Grad_Symbol(DependendSymbol):
4892        @type format: C{str}        @type format: C{str}
4893        @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.
4894        @rtype: C{str}        @rtype: C{str}
4895        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4896        """        """
4897        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4898           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4487  def integrate(arg,where=None): Line 4945  def integrate(arg,where=None):
4945                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4946      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4947      @return: integral of arg.      @return: integral of arg.
4948      @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}      @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4949      """      """
4950      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4951         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
# Line 4498  def integrate(arg,where=None): Line 4956  def integrate(arg,where=None):
4956         else:         else:
4957            return arg._integrate()            return arg._integrate()
4958      else:      else:
4959        raise TypeError,"integrate: Unknown argument type."         arg2=escript.Data(arg,where)
4960           if arg2.getRank()==0:
4961              return arg2._integrate()[0]
4962           else:
4963              return arg2._integrate()
4964    
4965  class Integrate_Symbol(DependendSymbol):  class Integrate_Symbol(DependendSymbol):
4966     """     """
# Line 4525  class Integrate_Symbol(DependendSymbol): Line 4987  class Integrate_Symbol(DependendSymbol):
4987        @type format: C{str}        @type format: C{str}
4988        @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.
4989        @rtype: C{str}        @rtype: C{str}
4990        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4991        """        """
4992        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4993           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4570  class Integrate_Symbol(DependendSymbol): Line 5032  class Integrate_Symbol(DependendSymbol):
5032    
5033  def interpolate(arg,where):  def interpolate(arg,where):
5034      """      """
5035      interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where. If the argument C{arg} has the requested function space
5036        C{where} no interpolation is performed and C{arg} is returned.
5037    
5038      @param arg: interpolant      @param arg: interpolant
5039      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
5040      @param where: FunctionSpace to be interpolated to      @param where: FunctionSpace to be interpolated to
5041      @type where: L{escript.FunctionSpace}      @type where: L{escript.FunctionSpace}
5042      @return: interpolated argument      @return: interpolated argument
5043      @rtype:  C{escript.Data} or L{Symbol}      @rtype: C{escript.Data} or L{Symbol}
5044      """      """
5045      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5046         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
5047      else:      else:
5048         return escript.Data(arg,where)         if where == arg.getFunctionSpace():
5049              return arg
5050           else:
5051              return escript.Data(arg,where)
5052    
5053  class Interpolate_Symbol(DependendSymbol):  class Interpolate_Symbol(DependendSymbol):
5054     """     """
# Line 4608  class Interpolate_Symbol(DependendSymbol Line 5074  class Interpolate_Symbol(DependendSymbol
5074        @type format: C{str}        @type format: C{str}
5075        @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.
5076        @rtype: C{str}        @rtype: C{str}
5077        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
5078        """        """
5079        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
5080           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4661  def div(arg,where=None): Line 5127  def div(arg,where=None):
5127                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
5128      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
5129      @return: divergence of arg.      @return: divergence of arg.
5130      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5131      """      """
5132      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5133          dim=arg.getDim()          dim=arg.getDim()
# Line 4683  def jump(arg,domain=None): Line 5149  def jump(arg,domain=None):
5149                     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.
5150      @type domain: C{None} or L{escript.Domain}      @type domain: C{None} or L{escript.Domain}
5151      @return: jump of arg      @return: jump of arg
5152      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5153      """      """
5154      if domain==None: domain=arg.getDomain()      if domain==None: domain=arg.getDomain()
5155      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
# Line 4695  def L2(arg): Line 5161  def L2(arg):
5161      @param arg: function which L2 to be calculated.      @param arg: function which L2 to be calculated.
5162      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
5163      @return: L2 norm of arg.      @return: L2 norm of arg.
5164      @rtype:  L{float} or L{Symbol}      @rtype: L{float} or L{Symbol}
5165      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5166      """      """
5167      return sqrt(integrate(inner(arg,arg)))      return sqrt(integrate(inner(arg,arg)))

Legend:
Removed from v.587  
changed lines
  Added in v.1357

  ViewVC Help
Powered by ViewVC 1.1.26