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

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

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

revision 429 by gross, Wed Jan 11 05:53:40 2006 UTC revision 795 by ksteube, Mon Jul 31 01:23:58 2006 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2    
3  """  """
4  Utility functions for escript  Utility functions for escript
5    
 @remark:  This module is under construction and is still tested!!!  
   
6  @var __author__: name of author  @var __author__: name of author
7  @var __licence__: licence agreement  @var __copyright__: copyrights
8    @var __license__: licence agreement
9  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
10  @var __version__: version  @var __version__: version
11  @var __date__: date of the version  @var __date__: date of the version
12  """  """
13                                                                                                                                                                                                                                                                                                                                                                                                            
14  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
15  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                        http://www.access.edu.au
17                    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19                 http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
21  __version__="$Revision$"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
# Line 33  import numarray Line 27  import numarray
27  import escript  import escript
28  import os  import os
29    
 # missing tests:  
   
 # def pokeShape(arg):  
 # def pokeDim(arg):  
 # def commonShape(arg0,arg1):  
 # def commonDim(*args):  
 # def testForZero(arg):  
 # def matchType(arg0=0.,arg1=0.):  
 # def matchShape(arg0,arg1):  
   
 # def transpose(arg,axis=None):  
 # def trace(arg,axis0=0,axis1=1):  
 # def reorderComponents(arg,index):  
   
 # def integrate(arg,where=None):  
 # def interpolate(arg,where):  
 # def div(arg,where=None):  
 # def grad(arg,where=None):  
   
 #  
 # slicing: get  
 #          set  
 #  
 # and derivatives  
   
30  #=========================================================  #=========================================================
31  #   some helpers:  #   some helpers:
32  #=========================================================  #=========================================================
# Line 65  def saveVTK(filename,domain=None,**data) Line 34  def saveVTK(filename,domain=None,**data)
34      """      """
35      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.
36    
37      Example:      Example::
38    
39         tmp=Scalar(..)         tmp=Scalar(..)
40         v=Vector(..)         v=Vector(..)
# Line 93  def saveDX(filename,domain=None,**data): Line 62  def saveDX(filename,domain=None,**data):
62      """      """
63      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.
64    
65      Example:      Example::
66    
67         tmp=Scalar(..)         tmp=Scalar(..)
68         v=Vector(..)         v=Vector(..)
# Line 122  def kronecker(d=3): Line 91  def kronecker(d=3):
91     return the kronecker S{delta}-symbol     return the kronecker S{delta}-symbol
92    
93     @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
94     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
95     @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
96     @rtype d: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2.
    @remark: the function is identical L{identity}  
97     """     """
98     return identityTensor(d)     return identityTensor(d)
99    
# Line 140  def identity(shape=()): Line 108  def identity(shape=()):
108     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
109     """     """
110     if len(shape)>0:     if len(shape)>0:
111        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
112        if len(shape)==1:        if len(shape)==1:
113            for i0 in range(shape[0]):            for i0 in range(shape[0]):
114               out[i0,i0]=1.               out[i0,i0]=1.
   
115        elif len(shape)==2:        elif len(shape)==2:
116            for i0 in range(shape[0]):            for i0 in range(shape[0]):
117               for i1 in range(shape[1]):               for i1 in range(shape[1]):
# Line 160  def identityTensor(d=3): Line 127  def identityTensor(d=3):
127     return the dxd identity matrix     return the dxd identity matrix
128    
129     @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
130     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
131     @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
132     @rtype: L{numarray.NumArray} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
133     """     """
134     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
135        d=d.getDim()         return escript.Data(identity((d.getDim(),)),d)
136     return identity(shape=(d,))     elif isinstance(d,escript.Domain):
137           return identity((d.getDim(),))
138       else:
139           return identity((d,))
140    
141  def identityTensor4(d=3):  def identityTensor4(d=3):
142     """     """
# Line 175  def identityTensor4(d=3): Line 145  def identityTensor4(d=3):
145     @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
146     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
147     @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
148     @rtype: L{numarray.NumArray} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
149     """     """
150     if hasattr(d,"getDim"):     if isinstance(d,escript.FunctionSpace):
151        d=d.getDim()         return escript.Data(identity((d.getDim(),d.getDim())),d)
152     return identity((d,d))     elif isinstance(d,escript.Domain):
153           return identity((d.getDim(),d.getDim()))
154       else:
155           return identity((d,d))
156    
157  def unitVector(i=0,d=3):  def unitVector(i=0,d=3):
158     """     """
# Line 188  def unitVector(i=0,d=3): Line 161  def unitVector(i=0,d=3):
161     @param i: index     @param i: index
162     @type i: C{int}     @type i: C{int}
163     @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
164     @type d: C{int} or any object with a C{getDim} method     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
165     @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
166     @rtype: L{numarray.NumArray} of rank 1.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
167     """     """
168     return kronecker(d)[i]     return kronecker(d)[i]
169    
# Line 246  def inf(arg): Line 219  def inf(arg):
219    
220      @param arg: argument      @param arg: argument
221      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
222      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
223      @rtype: C{float}      @rtype: C{float}
224      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
225      """      """
# Line 335  def commonDim(*args): Line 308  def commonDim(*args):
308      """      """
309      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.
310    
311      @param *args: given objects      @param args: given objects
312      @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
313               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
314      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 357  def testForZero(arg): Line 330  def testForZero(arg):
330    
331      @param arg: a given object      @param arg: a given object
332      @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}
333      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
334      @rtype : C{bool}      @rtype: C{bool}
335      """      """
336      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
337         return not Lsup(arg)>0.         return not Lsup(arg)>0.
# Line 391  def matchType(arg0=0.,arg1=0.): Line 364  def matchType(arg0=0.,arg1=0.):
364         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
365            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
366         elif isinstance(arg1,float):         elif isinstance(arg1,float):
367            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
368         elif isinstance(arg1,int):         elif isinstance(arg1,int):
369            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
370         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
371            pass            pass
372         else:         else:
# Line 417  def matchType(arg0=0.,arg1=0.): Line 390  def matchType(arg0=0.,arg1=0.):
390         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
391            pass            pass
392         elif isinstance(arg1,float):         elif isinstance(arg1,float):
393            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
394         elif isinstance(arg1,int):         elif isinstance(arg1,int):
395            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
396         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
397            pass            pass
398         else:         else:
399            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
400      elif isinstance(arg0,float):      elif isinstance(arg0,float):
401         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
402            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
403         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
404            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
405         elif isinstance(arg1,float):         elif isinstance(arg1,float):
406            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
407            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
408         elif isinstance(arg1,int):         elif isinstance(arg1,int):
409            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
410            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
411         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
412            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
413         else:         else:
414            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
415      elif isinstance(arg0,int):      elif isinstance(arg0,int):
416         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
417            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
418         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
419            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
420         elif isinstance(arg1,float):         elif isinstance(arg1,float):
421            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
422            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
423         elif isinstance(arg1,int):         elif isinstance(arg1,int):
424            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
425            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
426         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
427            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
428         else:         else:
429            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
430      else:      else:
# Line 461  def matchType(arg0=0.,arg1=0.): Line 434  def matchType(arg0=0.,arg1=0.):
434    
435  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
436      """      """
437            return representations of arg0 amd arg1 which ahve the same shape
438    
439      If shape is not given the shape "largest" shape of args is used.      @param arg0: a given object
440        @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
441      @param args: a given ob      @param arg1: a given object
442      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
443      @return: True if the argument is identical to zero.      @return: arg0 and arg1 where copies are returned when the shape has to be changed.
444      @rtype: C{list} of C{int}      @rtype: C{tuple}
445      """      """
446      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
447      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
448      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
449      if len(sh0)<len(sh):      if len(sh0)<len(sh):
450         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
451      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
452         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
453      else:      else:
454         return arg0,arg1         return arg0,arg1
455  #=========================================================  #=========================================================
# Line 496  class Symbol(object): Line 469  class Symbol(object):
469         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
470         symbols or any other object.         symbols or any other object.
471    
472         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
473         @type arg: C{list}         @type args: C{list}
474         @param shape: the shape         @param shape: the shape
475         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
476         @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
# Line 540  class Symbol(object): Line 513  class Symbol(object):
513         """         """
514         the shape of the symbol.         the shape of the symbol.
515    
516         @return : the shape of the symbol.         @return: the shape of the symbol.
517         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
518         """         """
519         return self.__shape         return self.__shape
# Line 549  class Symbol(object): Line 522  class Symbol(object):
522         """         """
523         the spatial dimension         the spatial dimension
524    
525         @return : the spatial dimension         @return: the spatial dimension
526         @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.
527         """         """
528         return self.__dim         return self.__dim
# Line 573  class Symbol(object): Line 546  class Symbol(object):
546         """         """
547         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.
548    
549         @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].
550         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
551         @rtype: C{list} of objects         @rtype: C{list} of objects
552         @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.
553         """         """
554         out=[]         out=[]
555         for a in self.getArgument():         for a in self.getArgument():
# Line 602  class Symbol(object): Line 575  class Symbol(object):
575            else:            else:
576                s=pokeShape(s)+arg.getShape()                s=pokeShape(s)+arg.getShape()
577                if len(s)>0:                if len(s)>0:
578                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
579                else:                else:
580                   out.append(a)                   out.append(a)
581         return out         return out
# Line 692  class Symbol(object): Line 665  class Symbol(object):
665         else:         else:
666            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
667            if len(s)>0:            if len(s)>0:
668               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
669            else:            else:
670               return 0.               return 0.
671    
# Line 700  class Symbol(object): Line 673  class Symbol(object):
673         """         """
674         returns -self.         returns -self.
675    
676         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
677         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
678         """         """
679         return self*(-1.)         return self*(-1.)
# Line 709  class Symbol(object): Line 682  class Symbol(object):
682         """         """
683         returns +self.         returns +self.
684    
685         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
686         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
687         """         """
688         return self*(1.)         return self*(1.)
689    
690     def __abs__(self):     def __abs__(self):
691         """         """
692         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
693         """         """
694         return Abs_Symbol(self)         return Abs_Symbol(self)
695    
# Line 726  class Symbol(object): Line 699  class Symbol(object):
699    
700         @param other: object to be added to this object         @param other: object to be added to this object
701         @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}.
702         @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}
703         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
704         """         """
705         return add(self,other)         return add(self,other)
# Line 737  class Symbol(object): Line 710  class Symbol(object):
710    
711         @param other: object this object is added to         @param other: object this object is added to
712         @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}.
713         @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
714         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
715         """         """
716         return add(other,self)         return add(other,self)
# Line 748  class Symbol(object): Line 721  class Symbol(object):
721    
722         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
723         @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}.
724         @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
725         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
726         """         """
727         return add(self,-other)         return add(self,-other)
# Line 759  class Symbol(object): Line 732  class Symbol(object):
732    
733         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
734         @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}.
735         @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}.
736         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
737         """         """
738         return add(-self,other)         return add(-self,other)
# Line 770  class Symbol(object): Line 743  class Symbol(object):
743    
744         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
745         @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}.
746         @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}.
747         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
748         """         """
749         return mult(self,other)         return mult(self,other)
# Line 781  class Symbol(object): Line 754  class Symbol(object):
754    
755         @param other: object this object is multiplied with         @param other: object this object is multiplied with
756         @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}.
757         @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.
758         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
759         """         """
760         return mult(other,self)         return mult(other,self)
# Line 792  class Symbol(object): Line 765  class Symbol(object):
765    
766         @param other: object dividing this object         @param other: object dividing this object
767         @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}.
768         @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}
769         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
770         """         """
771         return quotient(self,other)         return quotient(self,other)
# Line 803  class Symbol(object): Line 776  class Symbol(object):
776    
777         @param other: object dividing this object         @param other: object dividing this object
778         @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}.
779         @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
780         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
781         """         """
782         return quotient(other,self)         return quotient(other,self)
# Line 814  class Symbol(object): Line 787  class Symbol(object):
787    
788         @param other: exponent         @param other: exponent
789         @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}.
790         @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}
791         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
792         """         """
793         return power(self,other)         return power(self,other)
# Line 825  class Symbol(object): Line 798  class Symbol(object):
798    
799         @param other: basis         @param other: basis
800         @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}.
801         @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
802         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
803         """         """
804         return power(other,self)         return power(other,self)
805    
806       def __getitem__(self,index):
807           """
808           returns the slice defined by index
809    
810           @param index: defines a
811           @type index: C{slice} or C{int} or a C{tuple} of them
812           @return: a L{Symbol} representing the slice defined by index
813           @rtype: L{DependendSymbol}
814           """
815           return GetSlice_Symbol(self,index)
816    
817  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
818     """     """
819     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.
820     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  
821        
822     Example:     Example::
823        
824     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
825     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
826     print u1==u2       print u1==u2
827     False       False
828        
829        but       but::
830    
831     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
832     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
833     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
834     print u1==u2, u1==u3       print u1==u2, u1==u3
835     True False       True False
836    
837     @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.
838     """     """
# Line 880  class DependendSymbol(Symbol): Line 864  class DependendSymbol(Symbol):
864  #=========================================================  #=========================================================
865  #  Unary operations prserving the shape  #  Unary operations prserving the shape
866  #========================================================  #========================================================
867    class GetSlice_Symbol(DependendSymbol):
868       """
869       L{Symbol} representing getting a slice for a L{Symbol}
870       """
871       def __init__(self,arg,index):
872          """
873          initialization of wherePositive L{Symbol} with argument arg
874          @param arg: argument
875          @type arg: L{Symbol}.
876          @param index: defines index
877          @type index: C{slice} or C{int} or a C{tuple} of them
878          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
879          @raises ValueError: if a step is given
880          """
881          if not isinstance(index,tuple): index=(index,)
882          if len(index)>arg.getRank():
883               raise IndexError,"GetSlice_Symbol: index out of range."
884          sh=()
885          index2=()
886          for i in range(len(index)):
887             ix=index[i]
888             if isinstance(ix,int):
889                if ix<0 or ix>=arg.getShape()[i]:
890                   raise ValueError,"GetSlice_Symbol: index out of range."
891                index2=index2+(ix,)
892             else:
893               if not ix.step==None:
894                 raise ValueError,"GetSlice_Symbol: steping is not supported."
895               if ix.start==None:
896                  s=0
897               else:
898                  s=ix.start
899               if ix.stop==None:
900                  e=arg.getShape()[i]
901               else:
902                  e=ix.stop
903                  if e>arg.getShape()[i]:
904                     raise IndexError,"GetSlice_Symbol: index out of range."
905               index2=index2+(slice(s,e),)
906               if e>s:
907                   sh=sh+(e-s,)
908               elif s>e:
909                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
910          for i in range(len(index),arg.getRank()):
911              index2=index2+(slice(0,arg.getShape()[i]),)
912              sh=sh+(arg.getShape()[i],)
913          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
914    
915       def getMyCode(self,argstrs,format="escript"):
916          """
917          returns a program code that can be used to evaluate the symbol.
918    
919          @param argstrs: gives for each argument a string representing the argument for the evaluation.
920          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
921          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
922          @type format: C{str}
923          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
924          @rtype: C{str}
925          @raise NotImplementedError: if the requested format is not available
926          """
927          if format=="escript" or format=="str"  or format=="text":
928             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
929          else:
930             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
931    
932       def substitute(self,argvals):
933          """
934          assigns new values to symbols in the definition of the symbol.
935          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
936    
937          @param argvals: new values assigned to symbols
938          @type argvals: C{dict} with keywords of type L{Symbol}.
939          @return: result of the substitution process. Operations are executed as much as possible.
940          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
941          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
942          """
943          if argvals.has_key(self):
944             arg=argvals[self]
945             if self.isAppropriateValue(arg):
946                return arg
947             else:
948                raise TypeError,"%s: new value is not appropriate."%str(self)
949          else:
950             args=self.getSubstitutedArguments(argvals)
951             arg=args[0]
952             index=args[1]
953             return arg.__getitem__(index)
954    
955  def log10(arg):  def log10(arg):
956     """     """
957     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
958    
959     @param arg: argument     @param arg: argument
960     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
961     @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.
962     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
963     """     """
964     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 908  def wherePositive(arg): Line 980  def wherePositive(arg):
980    
981     @param arg: argument     @param arg: argument
982     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
983     @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.
984     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
985     """     """
986     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
987        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
988        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
989        return out        return out
990     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
991        return arg._wherePositive()        return arg._wherePositive()
# Line 954  class WherePositive_Symbol(DependendSymb Line 1026  class WherePositive_Symbol(DependendSymb
1026        @type format: C{str}        @type format: C{str}
1027        @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.
1028        @rtype: C{str}        @rtype: C{str}
1029        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1030        """        """
1031        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1032            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 990  def whereNegative(arg): Line 1062  def whereNegative(arg):
1062    
1063     @param arg: argument     @param arg: argument
1064     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1065     @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.
1066     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1067     """     """
1068     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1069        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1070        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1071        return out        return out
1072     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1073        return arg._whereNegative()        return arg._whereNegative()
# Line 1036  class WhereNegative_Symbol(DependendSymb Line 1108  class WhereNegative_Symbol(DependendSymb
1108        @type format: C{str}        @type format: C{str}
1109        @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.
1110        @rtype: C{str}        @rtype: C{str}
1111        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1112        """        """
1113        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1114            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1072  def whereNonNegative(arg): Line 1144  def whereNonNegative(arg):
1144    
1145     @param arg: argument     @param arg: argument
1146     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1147     @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.
1148     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1149     """     """
1150     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1151        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1152        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1153        return out        return out
1154     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1155        return arg._whereNonNegative()        return arg._whereNonNegative()
# Line 1102  def whereNonPositive(arg): Line 1174  def whereNonPositive(arg):
1174    
1175     @param arg: argument     @param arg: argument
1176     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1177     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1178     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1179     """     """
1180     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1181        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1182        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1183        return out        return out
1184     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1185        return arg._whereNonPositive()        return arg._whereNonPositive()
# Line 1134  def whereZero(arg,tol=0.): Line 1206  def whereZero(arg,tol=0.):
1206     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1207     @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.
1208     @type tol: C{float}     @type tol: C{float}
1209     @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.
1210     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1211     """     """
1212     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1213        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1214        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1215        return out        return out
1216     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1217        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1218     elif isinstance(arg,float):     elif isinstance(arg,float):
1219        if abs(arg)<=tol:        if abs(arg)<=tol:
1220          return 1.          return 1.
# Line 1183  class WhereZero_Symbol(DependendSymbol): Line 1252  class WhereZero_Symbol(DependendSymbol):
1252        @type format: C{str}        @type format: C{str}
1253        @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.
1254        @rtype: C{str}        @rtype: C{str}
1255        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1256        """        """
1257        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1258           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1217  def whereNonZero(arg,tol=0.): Line 1286  def whereNonZero(arg,tol=0.):
1286    
1287     @param arg: argument     @param arg: argument
1288     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1289     @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1290     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1291     """     """
1292     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1293        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1294        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1295        return out        return out
1296     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1297        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1298     elif isinstance(arg,float):     elif isinstance(arg,float):
1299        if abs(arg)>tol:        if abs(arg)>tol:
1300          return 1.          return 1.
# Line 1250  def sin(arg): Line 1316  def sin(arg):
1316    
1317     @param arg: argument     @param arg: argument
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     @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.
1320     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1321     """     """
1322     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1288  class Sin_Symbol(DependendSymbol): Line 1354  class Sin_Symbol(DependendSymbol):
1354        @type format: C{str}        @type format: C{str}
1355        @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.
1356        @rtype: C{str}        @rtype: C{str}
1357        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1358        """        """
1359        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1360            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1340  def cos(arg): Line 1406  def cos(arg):
1406    
1407     @param arg: argument     @param arg: argument
1408     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1409     @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.
1410     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1411     """     """
1412     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1378  class Cos_Symbol(DependendSymbol): Line 1444  class Cos_Symbol(DependendSymbol):
1444        @type format: C{str}        @type format: C{str}
1445        @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.
1446        @rtype: C{str}        @rtype: C{str}
1447        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1448        """        """
1449        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1450            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1430  def tan(arg): Line 1496  def tan(arg):
1496    
1497     @param arg: argument     @param arg: argument
1498     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1499     @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.
1500     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1501     """     """
1502     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1468  class Tan_Symbol(DependendSymbol): Line 1534  class Tan_Symbol(DependendSymbol):
1534        @type format: C{str}        @type format: C{str}
1535        @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.
1536        @rtype: C{str}        @rtype: C{str}
1537        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1538        """        """
1539        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1540            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1520  def asin(arg): Line 1586  def asin(arg):
1586    
1587     @param arg: argument     @param arg: argument
1588     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1589     @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.
1590     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1591     """     """
1592     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1558  class Asin_Symbol(DependendSymbol): Line 1624  class Asin_Symbol(DependendSymbol):
1624        @type format: C{str}        @type format: C{str}
1625        @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.
1626        @rtype: C{str}        @rtype: C{str}
1627        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1628        """        """
1629        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1630            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1610  def acos(arg): Line 1676  def acos(arg):
1676    
1677     @param arg: argument     @param arg: argument
1678     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1679     @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.
1680     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1681     """     """
1682     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1648  class Acos_Symbol(DependendSymbol): Line 1714  class Acos_Symbol(DependendSymbol):
1714        @type format: C{str}        @type format: C{str}
1715        @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.
1716        @rtype: C{str}        @rtype: C{str}
1717        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1718        """        """
1719        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1720            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1700  def atan(arg): Line 1766  def atan(arg):
1766    
1767     @param arg: argument     @param arg: argument
1768     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1769     @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.
1770     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1771     """     """
1772     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1738  class Atan_Symbol(DependendSymbol): Line 1804  class Atan_Symbol(DependendSymbol):
1804        @type format: C{str}        @type format: C{str}
1805        @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.
1806        @rtype: C{str}        @rtype: C{str}
1807        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1808        """        """
1809        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1810            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1790  def sinh(arg): Line 1856  def sinh(arg):
1856    
1857     @param arg: argument     @param arg: argument
1858     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1859     @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.
1860     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1861     """     """
1862     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1828  class Sinh_Symbol(DependendSymbol): Line 1894  class Sinh_Symbol(DependendSymbol):
1894        @type format: C{str}        @type format: C{str}
1895        @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.
1896        @rtype: C{str}        @rtype: C{str}
1897        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1898        """        """
1899        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1900            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1880  def cosh(arg): Line 1946  def cosh(arg):
1946    
1947     @param arg: argument     @param arg: argument
1948     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1949     @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.
1950     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1951     """     """
1952     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1918  class Cosh_Symbol(DependendSymbol): Line 1984  class Cosh_Symbol(DependendSymbol):
1984        @type format: C{str}        @type format: C{str}
1985        @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.
1986        @rtype: C{str}        @rtype: C{str}
1987        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1988        """        """
1989        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1990            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1970  def tanh(arg): Line 2036  def tanh(arg):
2036    
2037     @param arg: argument     @param arg: argument
2038     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2039     @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.
2040     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2041     """     """
2042     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2008  class Tanh_Symbol(DependendSymbol): Line 2074  class Tanh_Symbol(DependendSymbol):
2074        @type format: C{str}        @type format: C{str}
2075        @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.
2076        @rtype: C{str}        @rtype: C{str}
2077        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2078        """        """
2079        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2080            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2060  def asinh(arg): Line 2126  def asinh(arg):
2126    
2127     @param arg: argument     @param arg: argument
2128     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2129     @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.
2130     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2131     """     """
2132     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2098  class Asinh_Symbol(DependendSymbol): Line 2164  class Asinh_Symbol(DependendSymbol):
2164        @type format: C{str}        @type format: C{str}
2165        @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.
2166        @rtype: C{str}        @rtype: C{str}
2167        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2168        """        """
2169        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2170            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2150  def acosh(arg): Line 2216  def acosh(arg):
2216    
2217     @param arg: argument     @param arg: argument
2218     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2219     @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.
2220     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2221     """     """
2222     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2188  class Acosh_Symbol(DependendSymbol): Line 2254  class Acosh_Symbol(DependendSymbol):
2254        @type format: C{str}        @type format: C{str}
2255        @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.
2256        @rtype: C{str}        @rtype: C{str}
2257        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2258        """        """
2259        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2260            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2240  def atanh(arg): Line 2306  def atanh(arg):
2306    
2307     @param arg: argument     @param arg: argument
2308     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2309     @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.
2310     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2311     """     """
2312     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2278  class Atanh_Symbol(DependendSymbol): Line 2344  class Atanh_Symbol(DependendSymbol):
2344        @type format: C{str}        @type format: C{str}
2345        @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.
2346        @rtype: C{str}        @rtype: C{str}
2347        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2348        """        """
2349        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2350            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2330  def exp(arg): Line 2396  def exp(arg):
2396    
2397     @param arg: argument     @param arg: argument
2398     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2399     @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.
2400     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2401     """     """
2402     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2368  class Exp_Symbol(DependendSymbol): Line 2434  class Exp_Symbol(DependendSymbol):
2434        @type format: C{str}        @type format: C{str}
2435        @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.
2436        @rtype: C{str}        @rtype: C{str}
2437        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2438        """        """
2439        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2440            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2420  def sqrt(arg): Line 2486  def sqrt(arg):
2486    
2487     @param arg: argument     @param arg: argument
2488     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2489     @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.
2490     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2491     """     """
2492     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2458  class Sqrt_Symbol(DependendSymbol): Line 2524  class Sqrt_Symbol(DependendSymbol):
2524        @type format: C{str}        @type format: C{str}
2525        @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.
2526        @rtype: C{str}        @rtype: C{str}
2527        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2528        """        """
2529        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2530            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2510  def log(arg): Line 2576  def log(arg):
2576    
2577     @param arg: argument     @param arg: argument
2578     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2579     @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.
2580     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2581     """     """
2582     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2548  class Log_Symbol(DependendSymbol): Line 2614  class Log_Symbol(DependendSymbol):
2614        @type format: C{str}        @type format: C{str}
2615        @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.
2616        @rtype: C{str}        @rtype: C{str}
2617        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2618        """        """
2619        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2620            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2600  def sign(arg): Line 2666  def sign(arg):
2666    
2667     @param arg: argument     @param arg: argument
2668     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2669     @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.
2670     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2671     """     """
2672     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2648  class Abs_Symbol(DependendSymbol): Line 2714  class Abs_Symbol(DependendSymbol):
2714        @type format: C{str}        @type format: C{str}
2715        @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.
2716        @rtype: C{str}        @rtype: C{str}
2717        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2718        """        """
2719        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2720            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2700  def minval(arg): Line 2766  def minval(arg):
2766    
2767     @param arg: argument     @param arg: argument
2768     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2769     @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.
2770     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2771     """     """
2772     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2741  class Minval_Symbol(DependendSymbol): Line 2807  class Minval_Symbol(DependendSymbol):
2807        @type format: C{str}        @type format: C{str}
2808        @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.
2809        @rtype: C{str}        @rtype: C{str}
2810        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2811        """        """
2812        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2813            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2777  def maxval(arg): Line 2843  def maxval(arg):
2843    
2844     @param arg: argument     @param arg: argument
2845     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2846     @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.
2847     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2848     """     """
2849     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2818  class Maxval_Symbol(DependendSymbol): Line 2884  class Maxval_Symbol(DependendSymbol):
2884        @type format: C{str}        @type format: C{str}
2885        @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.
2886        @rtype: C{str}        @rtype: C{str}
2887        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2888        """        """
2889        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2890            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2854  def length(arg): Line 2920  def length(arg):
2920    
2921     @param arg: argument     @param arg: argument
2922     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2923     @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.
2924     """     """
2925     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
2926    
# Line 2873  def trace(arg,axis_offset=0): Line 2939  def trace(arg,axis_offset=0):
2939     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
2940        sh=arg.shape        sh=arg.shape
2941        if len(sh)<2:        if len(sh)<2:
2942          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2943        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
2944          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
2945        s1=1        s1=1
2946        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
2947        s2=1        s2=1
2948        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2949        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
2950          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)
2951        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))
2952        out=numarray.zeros([s1,s2],numarray.Float)        out=numarray.zeros([s1,s2],numarray.Float64)
2953        for i1 in range(s1):        for i1 in range(s1):
2954          for i2 in range(s2):          for i2 in range(s2):
2955              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
2956        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2957        return out        return out
2958     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
2959        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
2960            raise ValueError,"rank of argument must be greater than 1"
2961          if axis_offset<0 or axis_offset>arg.getRank()-2:
2962            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
2963          s=list(arg.getShape())        
2964          if not s[axis_offset] == s[axis_offset+1]:
2965            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2966          return arg._matrixtrace(axis_offset)
2967     elif isinstance(arg,float):     elif isinstance(arg,float):
2968        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
2969     elif isinstance(arg,int):     elif isinstance(arg,int):
2970        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
2971     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
2972        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
2973     else:     else:
2974        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
2975    
 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  
2976  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
2977     """     """
2978     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 2953  class Trace_Symbol(DependendSymbol): Line 2987  class Trace_Symbol(DependendSymbol):
2987        @type axis_offset: C{int}        @type axis_offset: C{int}
2988        """        """
2989        if arg.getRank()<2:        if arg.getRank()<2:
2990          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2991        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
2992          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
2993        s=list(arg.getShape())                s=list(arg.getShape())        
2994        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
2995          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)
2996        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())
2997    
2998     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 2971  class Trace_Symbol(DependendSymbol): Line 3005  class Trace_Symbol(DependendSymbol):
3005        @type format: C{str}        @type format: C{str}
3006        @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.
3007        @rtype: C{str}        @rtype: C{str}
3008        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3009        """        """
3010        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3011           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3013  class Trace_Symbol(DependendSymbol): Line 3047  class Trace_Symbol(DependendSymbol):
3047        else:        else:
3048           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3049    
3050    def transpose(arg,axis_offset=None):
3051       """
3052       returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.
3053    
3054       @param arg: argument
3055       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3056       @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.
3057                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3058       @type axis_offset: C{int}
3059       @return: transpose of arg
3060       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3061       """
3062       if isinstance(arg,numarray.NumArray):
3063          if axis_offset==None: axis_offset=int(arg.rank/2)
3064          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3065       elif isinstance(arg,escript.Data):
3066          r=arg.getRank()
3067          if axis_offset==None: axis_offset=int(r/2)
3068          if axis_offset<0 or axis_offset>r:
3069            raise ValueError,"axis_offset must be between 0 and %s"%r
3070          return arg._transpose(axis_offset)
3071       elif isinstance(arg,float):
3072          if not ( axis_offset==0 or axis_offset==None):
3073            raise ValueError,"axis_offset must be 0 for float argument"
3074          return arg
3075       elif isinstance(arg,int):
3076          if not ( axis_offset==0 or axis_offset==None):
3077            raise ValueError,"axis_offset must be 0 for int argument"
3078          return float(arg)
3079       elif isinstance(arg,Symbol):
3080          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3081          return Transpose_Symbol(arg,axis_offset)
3082       else:
3083          raise TypeError,"Unknown argument type."
3084    
3085    class Transpose_Symbol(DependendSymbol):
3086       """
3087       L{Symbol} representing the result of the transpose function
3088       """
3089       def __init__(self,arg,axis_offset=None):
3090          """
3091          initialization of transpose L{Symbol} with argument arg
3092    
3093          @param arg: argument of function
3094          @type arg: L{Symbol}.
3095          @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.
3096                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3097          @type axis_offset: C{int}
3098          """
3099          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3100          if axis_offset<0 or axis_offset>arg.getRank():
3101            raise ValueError,"axis_offset must be between 0 and %s"%r
3102          s=arg.getShape()
3103          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3104    
3105       def getMyCode(self,argstrs,format="escript"):
3106          """
3107          returns a program code that can be used to evaluate the symbol.
3108    
3109          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3110          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3111          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3112          @type format: C{str}
3113          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3114          @rtype: C{str}
3115          @raise NotImplementedError: if the requested format is not available
3116          """
3117          if format=="escript" or format=="str"  or format=="text":
3118             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3119          else:
3120             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3121    
3122       def substitute(self,argvals):
3123          """
3124          assigns new values to symbols in the definition of the symbol.
3125          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3126    
3127          @param argvals: new values assigned to symbols
3128          @type argvals: C{dict} with keywords of type L{Symbol}.
3129          @return: result of the substitution process. Operations are executed as much as possible.
3130          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3131          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3132          """
3133          if argvals.has_key(self):
3134             arg=argvals[self]
3135             if self.isAppropriateValue(arg):
3136                return arg
3137             else:
3138                raise TypeError,"%s: new value is not appropriate."%str(self)
3139          else:
3140             arg=self.getSubstitutedArguments(argvals)
3141             return transpose(arg[0],axis_offset=arg[1])
3142    
3143       def diff(self,arg):
3144          """
3145          differential of this object
3146    
3147          @param arg: the derivative is calculated with respect to arg
3148          @type arg: L{escript.Symbol}
3149          @return: derivative with respect to C{arg}
3150          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3151          """
3152          if arg==self:
3153             return identity(self.getShape())
3154          else:
3155             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3156    def symmetric(arg):
3157        """
3158        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3159    
3160        @param arg: square matrix. Must have rank 2 or 4 and be square.
3161        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3162        @return: symmetric part of arg
3163        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3164        """
3165        if isinstance(arg,numarray.NumArray):
3166          if arg.rank==2:
3167            if not (arg.shape[0]==arg.shape[1]):
3168               raise ValueError,"argument must be square."
3169          elif arg.rank==4:
3170            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3171               raise ValueError,"argument must be square."
3172          else:
3173            raise ValueError,"rank 2 or 4 is required."
3174          return (arg+transpose(arg))/2
3175        elif isinstance(arg,escript.Data):
3176          if arg.getRank()==2:
3177            if not (arg.getShape()[0]==arg.getShape()[1]):
3178               raise ValueError,"argument must be square."
3179            return arg._symmetric()
3180          elif arg.getRank()==4:
3181            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3182               raise ValueError,"argument must be square."
3183            return arg._symmetric()
3184          else:
3185            raise ValueError,"rank 2 or 4 is required."
3186        elif isinstance(arg,float):
3187          return arg
3188        elif isinstance(arg,int):
3189          return float(arg)
3190        elif isinstance(arg,Symbol):
3191          if arg.getRank()==2:
3192            if not (arg.getShape()[0]==arg.getShape()[1]):
3193               raise ValueError,"argument must be square."
3194          elif arg.getRank()==4:
3195            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3196               raise ValueError,"argument must be square."
3197          else:
3198            raise ValueError,"rank 2 or 4 is required."
3199          return (arg+transpose(arg))/2
3200        else:
3201          raise TypeError,"symmetric: Unknown argument type."
3202    
3203    def nonsymmetric(arg):
3204        """
3205        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3206    
3207        @param arg: square matrix. Must have rank 2 or 4 and be square.
3208        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3209        @return: nonsymmetric part of arg
3210        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3211        """
3212        if isinstance(arg,numarray.NumArray):
3213          if arg.rank==2:
3214            if not (arg.shape[0]==arg.shape[1]):
3215               raise ValueError,"nonsymmetric: argument must be square."
3216          elif arg.rank==4:
3217            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3218               raise ValueError,"nonsymmetric: argument must be square."
3219          else:
3220            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3221          return (arg-transpose(arg))/2
3222        elif isinstance(arg,escript.Data):
3223          if arg.getRank()==2:
3224            if not (arg.getShape()[0]==arg.getShape()[1]):
3225               raise ValueError,"argument must be square."
3226            return arg._nonsymmetric()
3227          elif arg.getRank()==4:
3228            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3229               raise ValueError,"argument must be square."
3230            return arg._nonsymmetric()
3231          else:
3232            raise ValueError,"rank 2 or 4 is required."
3233        elif isinstance(arg,float):
3234          return arg
3235        elif isinstance(arg,int):
3236          return float(arg)
3237        elif isinstance(arg,Symbol):
3238          if arg.getRank()==2:
3239            if not (arg.getShape()[0]==arg.getShape()[1]):
3240               raise ValueError,"nonsymmetric: argument must be square."
3241          elif arg.getRank()==4:
3242            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3243               raise ValueError,"nonsymmetric: argument must be square."
3244          else:
3245            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3246          return (arg-transpose(arg))/2
3247        else:
3248          raise TypeError,"nonsymmetric: Unknown argument type."
3249    
3250    def inverse(arg):
3251        """
3252        returns the inverse of the square matrix arg.
3253    
3254        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3255        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3256        @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3257        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3258        @note: for L{escript.Data} objects the dimension is restricted to 3.
3259        """
3260        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3261        if isinstance(arg,numarray.NumArray):
3262          return numarray.linear_algebra.inverse(arg)
3263        elif isinstance(arg,escript.Data):
3264          return escript_inverse(arg)
3265        elif isinstance(arg,float):
3266          return 1./arg
3267        elif isinstance(arg,int):
3268          return 1./float(arg)
3269        elif isinstance(arg,Symbol):
3270          return Inverse_Symbol(arg)
3271        else:
3272          raise TypeError,"inverse: Unknown argument type."
3273    
3274    def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3275          "arg is a Data objects!!!"
3276          if not arg.getRank()==2:
3277            raise ValueError,"escript_inverse: argument must have rank 2"
3278          s=arg.getShape()      
3279          if not s[0] == s[1]:
3280            raise ValueError,"escript_inverse: argument must be a square matrix."
3281          out=escript.Data(0.,s,arg.getFunctionSpace())
3282          if s[0]==1:
3283              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3284                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3285              out[0,0]=1./arg[0,0]
3286          elif s[0]==2:
3287              A11=arg[0,0]
3288              A12=arg[0,1]
3289              A21=arg[1,0]
3290              A22=arg[1,1]
3291              D = A11*A22-A12*A21
3292              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3293                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3294              D=1./D
3295              out[0,0]= A22*D
3296              out[1,0]=-A21*D
3297              out[0,1]=-A12*D
3298              out[1,1]= A11*D
3299          elif s[0]==3:
3300              A11=arg[0,0]
3301              A21=arg[1,0]
3302              A31=arg[2,0]
3303              A12=arg[0,1]
3304              A22=arg[1,1]
3305              A32=arg[2,1]
3306              A13=arg[0,2]
3307              A23=arg[1,2]
3308              A33=arg[2,2]
3309              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3310              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3311                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3312              D=1./D
3313              out[0,0]=(A22*A33-A23*A32)*D
3314              out[1,0]=(A31*A23-A21*A33)*D
3315              out[2,0]=(A21*A32-A31*A22)*D
3316              out[0,1]=(A13*A32-A12*A33)*D
3317              out[1,1]=(A11*A33-A31*A13)*D
3318              out[2,1]=(A12*A31-A11*A32)*D
3319              out[0,2]=(A12*A23-A13*A22)*D
3320              out[1,2]=(A13*A21-A11*A23)*D
3321              out[2,2]=(A11*A22-A12*A21)*D
3322          else:
3323             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3324          return out
3325    
3326    class Inverse_Symbol(DependendSymbol):
3327       """
3328       L{Symbol} representing the result of the inverse function
3329       """
3330       def __init__(self,arg):
3331          """
3332          initialization of inverse L{Symbol} with argument arg
3333          @param arg: argument of function
3334          @type arg: L{Symbol}.
3335          """
3336          if not arg.getRank()==2:
3337            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3338          s=arg.getShape()
3339          if not s[0] == s[1]:
3340            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3341          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3342    
3343       def getMyCode(self,argstrs,format="escript"):
3344          """
3345          returns a program code that can be used to evaluate the symbol.
3346    
3347          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3348          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3349          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3350          @type format: C{str}
3351          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3352          @rtype: C{str}
3353          @raise NotImplementedError: if the requested format is not available
3354          """
3355          if format=="escript" or format=="str"  or format=="text":
3356             return "inverse(%s)"%argstrs[0]
3357          else:
3358             raise NotImplementedError,"Inverse_Symbol does not provide program code for format %s."%format
3359    
3360       def substitute(self,argvals):
3361          """
3362          assigns new values to symbols in the definition of the symbol.
3363          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3364    
3365          @param argvals: new values assigned to symbols
3366          @type argvals: C{dict} with keywords of type L{Symbol}.
3367          @return: result of the substitution process. Operations are executed as much as possible.
3368          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3369          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3370          """
3371          if argvals.has_key(self):
3372             arg=argvals[self]
3373             if self.isAppropriateValue(arg):
3374                return arg
3375             else:
3376                raise TypeError,"%s: new value is not appropriate."%str(self)
3377          else:
3378             arg=self.getSubstitutedArguments(argvals)
3379             return inverse(arg[0])
3380    
3381       def diff(self,arg):
3382          """
3383          differential of this object
3384    
3385          @param arg: the derivative is calculated with respect to arg
3386          @type arg: L{escript.Symbol}
3387          @return: derivative with respect to C{arg}
3388          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3389          """
3390          if arg==self:
3391             return identity(self.getShape())
3392          else:
3393             return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3394    
3395    def eigenvalues(arg):
3396        """
3397        returns the eigenvalues of the square matrix arg.
3398    
3399        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3400                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3401        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3402        @return: the eigenvalues in increasing order.
3403        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3404        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3405        """
3406        if isinstance(arg,numarray.NumArray):
3407          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3408          out.sort()
3409          return out
3410        elif isinstance(arg,escript.Data):
3411          return arg._eigenvalues()
3412        elif isinstance(arg,Symbol):
3413          if not arg.getRank()==2:
3414            raise ValueError,"eigenvalues: argument must have rank 2"
3415          s=arg.getShape()      
3416          if not s[0] == s[1]:
3417            raise ValueError,"eigenvalues: argument must be a square matrix."
3418          if s[0]==1:
3419              return arg[0]
3420          elif s[0]==2:
3421              arg1=symmetric(arg)
3422              A11=arg1[0,0]
3423              A12=arg1[0,1]
3424              A22=arg1[1,1]
3425              trA=(A11+A22)/2.
3426              A11-=trA
3427              A22-=trA
3428              s=sqrt(A12**2-A11*A22)
3429              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3430          elif s[0]==3:
3431              arg1=symmetric(arg)
3432              A11=arg1[0,0]
3433              A12=arg1[0,1]
3434              A22=arg1[1,1]
3435              A13=arg1[0,2]
3436              A23=arg1[1,2]
3437              A33=arg1[2,2]
3438              trA=(A11+A22+A33)/3.
3439              A11-=trA
3440              A22-=trA
3441              A33-=trA
3442              A13_2=A13**2
3443              A23_2=A23**2
3444              A12_2=A12**2
3445              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3446              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3447              sq_p=sqrt(p/3.)
3448              alpha_3=acos(clip(-q*(sq_p+whereZero(p,0.)*1.e-15)**(-3.)/2.,-1.,1.))/3.  # whereZero is protection against divison by zero
3449              sq_p*=2.
3450              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3451               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3452               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3453              return trA+sq_p*f
3454          else:
3455             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3456        elif isinstance(arg,float):
3457          return arg
3458        elif isinstance(arg,int):
3459          return float(arg)
3460        else:
3461          raise TypeError,"eigenvalues: Unknown argument type."
3462    
3463    def eigenvalues_and_eigenvectors(arg):
3464        """
3465        returns the eigenvalues and eigenvectors of the square matrix arg.
3466    
3467        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3468                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3469        @type arg: L{escript.Data}
3470        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3471                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3472                 the eigenvector coresponding to the i-th eigenvalue.
3473        @rtype: L{tuple} of L{escript.Data}.
3474        @note: The dimension is restricted to 3.
3475        """
3476        if isinstance(arg,numarray.NumArray):
3477          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3478        elif isinstance(arg,escript.Data):
3479          return arg._eigenvalues_and_eigenvectors()
3480        elif isinstance(arg,Symbol):
3481          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3482        elif isinstance(arg,float):
3483          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3484        elif isinstance(arg,int):
3485          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3486        else:
3487          raise TypeError,"eigenvalues: Unknown argument type."
3488  #=======================================================  #=======================================================
3489  #  Binary operations:  #  Binary operations:
3490  #=======================================================  #=======================================================
# Line 3072  class Add_Symbol(DependendSymbol): Line 3544  class Add_Symbol(DependendSymbol):
3544        @type format: C{str}        @type format: C{str}
3545        @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.
3546        @rtype: C{str}        @rtype: C{str}
3547        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3548        """        """
3549        if format=="str" or format=="text":        if format=="str" or format=="text":
3550           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 3131  def mult(arg0,arg1): Line 3603  def mult(arg0,arg1):
3603         """         """
3604         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3605         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3606            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3607         else:         else:
3608            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3609                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3171  class Mult_Symbol(DependendSymbol): Line 3643  class Mult_Symbol(DependendSymbol):
3643        @type format: C{str}        @type format: C{str}
3644        @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.
3645        @rtype: C{str}        @rtype: C{str}
3646        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3647        """        """
3648        if format=="str" or format=="text":        if format=="str" or format=="text":
3649           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3231  def quotient(arg0,arg1): Line 3703  def quotient(arg0,arg1):
3703         """         """
3704         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3705         if testForZero(args[0]):         if testForZero(args[0]):
3706            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3707         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3708            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3709               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3276  class Quotient_Symbol(DependendSymbol): Line 3748  class Quotient_Symbol(DependendSymbol):
3748        @type format: C{str}        @type format: C{str}
3749        @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.
3750        @rtype: C{str}        @rtype: C{str}
3751        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3752        """        """
3753        if format=="str" or format=="text":        if format=="str" or format=="text":
3754           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3337  def power(arg0,arg1): Line 3809  def power(arg0,arg1):
3809         """         """
3810         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3811         if testForZero(args[0]):         if testForZero(args[0]):
3812            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3813         elif testForZero(args[1]):         elif testForZero(args[1]):
3814            return numarray.ones(args[0],numarray.Float)            return numarray.ones(pokeShape(args[1]),numarray.Float64)
3815         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3816            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3817         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 3380  class Power_Symbol(DependendSymbol): Line 3852  class Power_Symbol(DependendSymbol):
3852        @type format: C{str}        @type format: C{str}
3853        @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.
3854        @rtype: C{str}        @rtype: C{str}
3855        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3856        """        """
3857        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
3858           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3469  def clip(arg,minval=0.,maxval=1.): Line 3941  def clip(arg,minval=0.,maxval=1.):
3941      @param arg: argument      @param arg: argument
3942      @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}
3943      @param minval: lower range      @param minval: lower range
3944      @type arg: C{float}      @type minval: C{float}
3945      @param maxval: upper range      @param maxval: upper range
3946      @type arg: C{float}      @type maxval: C{float}
3947      @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
3948               less then maxval are unchanged.               less then maxval are unchanged.
3949      @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
# Line 3496  def inner(arg0,arg1): Line 3968  def inner(arg0,arg1):
3968      @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}
3969      @param arg1: second argument      @param arg1: second argument
3970      @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}
3971      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
3972      @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
3973      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
3974      """      """
# Line 3506  def inner(arg0,arg1): Line 3978  def inner(arg0,arg1):
3978          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
3979      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
3980    
3981    def outer(arg0,arg1):
3982        """
3983        the outer product of the two argument:
3984    
3985        out[t,s]=arg0[t]*arg1[s]
3986    
3987        where
3988    
3989            - s runs through arg0.Shape
3990            - t runs through arg1.Shape
3991    
3992        @param arg0: first argument
3993        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
3994        @param arg1: second argument
3995        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
3996        @return: the outer product of arg0 and arg1 at each data point
3997        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3998        """
3999        return generalTensorProduct(arg0,arg1,axis_offset=0)
4000    
4001  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4002      """      """
4003        see L{matrix_mult}
4004        """
4005        return matrix_mult(arg0,arg1)
4006    
4007    def matrix_mult(arg0,arg1):
4008        """
4009      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4010    
4011      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4012    
4013            or      or
4014    
4015      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4016    
4017      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.
4018    
4019      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4020      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3534  def matrixmult(arg0,arg1): Line 4032  def matrixmult(arg0,arg1):
4032          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4033      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4034    
4035  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4036      """      """
4037      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  
4038      """      """
4039      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
   
4040    
4041  def tensormult(arg0,arg1):  def tensor_mult(arg0,arg1):
4042      """      """
4043      the tensor product of the two argument:      the tensor product of the two argument:
   
4044            
4045      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4046    
4047      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4048    
4049                   or      or
4050    
4051      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4052    
# Line 3571  def tensormult(arg0,arg1): Line 4055  def tensormult(arg0,arg1):
4055    
4056      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]
4057                                
4058                   or      or
4059    
4060      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]
4061    
4062                   or      or
4063    
4064      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]
4065    
4066      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  
4067      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.
4068    
4069      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4070      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3596  def tensormult(arg0,arg1): Line 4080  def tensormult(arg0,arg1):
4080      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):
4081         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4082      else:      else:
4083          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4084    
4085  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4086      """      """
# Line 3604  def generalTensorProduct(arg0,arg1,axis_ Line 4088  def generalTensorProduct(arg0,arg1,axis_
4088    
4089      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4090    
4091      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:]  
4092    
4093      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]
4094      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4095            - t runs through arg1.Shape[axis_offset:]
4096    
4097      @param arg0: first argument      @param arg0: first argument
4098      @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}
4099      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4100      @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}
4101      @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.
4102      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
# Line 3626  def generalTensorProduct(arg0,arg1,axis_ Line 4109  def generalTensorProduct(arg0,arg1,axis_
4109             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4110         else:         else:
4111             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4112                 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)
4113             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4114             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4115             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 3636  def generalTensorProduct(arg0,arg1,axis_ Line 4119  def generalTensorProduct(arg0,arg1,axis_
4119             for i in sh1[:axis_offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4120             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4121             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4122             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4123             for i0 in range(d0):             for i0 in range(d0):
4124                      for i1 in range(d1):                      for i1 in range(d1):
4125                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
# Line 3652  def generalTensorProduct(arg0,arg1,axis_ Line 4135  def generalTensorProduct(arg0,arg1,axis_
4135                                    
4136  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4137     """     """
4138     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4139     """     """
4140     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4141         """         """
4142         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4143    
4144         @param arg0: numerator         @param arg0: first argument
4145         @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}.
4146         @param arg1: denominator         @param arg1: second argument
4147         @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}.
4148         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4149         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4150         """         """
4151         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
# Line 3685  class GeneralTensorProduct_Symbol(Depend Line 4168  class GeneralTensorProduct_Symbol(Depend
4168        @type format: C{str}        @type format: C{str}
4169        @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.
4170        @rtype: C{str}        @rtype: C{str}
4171        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4172        """        """
4173        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4174           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
# Line 3713  class GeneralTensorProduct_Symbol(Depend Line 4196  class GeneralTensorProduct_Symbol(Depend
4196           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4197           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4198    
4199    def escript_generalTensorProductNew(arg0,arg1,axis_offset):
4200        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4201        # calculate the return shape:
4202        shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4203        shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4204        shape10=arg1.getShape()[:axis_offset]
4205        shape1=arg1.getShape()[axis_offset:]
4206        if not shape01==shape10:
4207            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)
4208        # Figure out which functionspace to use (look at where operator+ is defined maybe in BinaryOp.h to get the logic for this)
4209        # fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4210        out=GeneralTensorProduct(arg0, arg1, axis_offset)
4211        return out
4212    
4213  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct
4214      "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!!!"
4215      # calculate the return shape:      # calculate the return shape:
# Line 3756  def escript_generalTensorProduct(arg0,ar Line 4253  def escript_generalTensorProduct(arg0,ar
4253      return out      return out
4254    
4255    
4256    def transposed_matrix_mult(arg0,arg1):
4257        """
4258        transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4259    
4260        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
4261    
4262        or
4263    
4264        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4265    
4266        The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4267    
4268        The first dimension of arg0 and arg1 must match.
4269    
4270        @param arg0: first argument of rank 2
4271        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4272        @param arg1: second argument of at least rank 1
4273        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4274        @return: the product of the transposed of arg0 and arg1 at each data point
4275        @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
4277        """
4278        sh0=pokeShape(arg0)
4279        sh1=pokeShape(arg1)
4280        if not len(sh0)==2 :
4281            raise ValueError,"first argument must have rank 2"
4282        if not len(sh1)==2 and not len(sh1)==1:
4283            raise ValueError,"second argument must have rank 1 or 2"
4284        return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4285    
4286    def transposed_tensor_mult(arg0,arg1):
4287        """
4288        the tensor product of the transposed of the first and the second argument
4289        
4290        for arg0 of rank 2 this is
4291    
4292        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4293    
4294        or
4295    
4296        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4297    
4298      
4299        and for arg0 of rank 4 this is
4300    
4301        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4302                  
4303        or
4304    
4305        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4306    
4307        or
4308    
4309        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4310    
4311        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4312        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4313    
4314        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4315    
4316        @param arg0: first argument of rank 2 or 4
4317        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4318        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4319        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4320        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4321        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4322        """
4323        sh0=pokeShape(arg0)
4324        sh1=pokeShape(arg1)
4325        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4326           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4327        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4328           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4329        else:
4330            raise ValueError,"first argument must have rank 2 or 4"
4331    
4332    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4333        """
4334        generalized tensor product of transposed of arg0 and arg1:
4335    
4336        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4337    
4338        where
4339    
4340            - s runs through arg0.Shape[axis_offset:]
4341            - r runs trough arg0.Shape[:axis_offset]
4342            - t runs through arg1.Shape[axis_offset:]
4343    
4344        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4345        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4346    
4347        @param arg0: first argument
4348        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4349        @param arg1: second argument
4350        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4351        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4352        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4353        """
4354        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4355        arg0,arg1=matchType(arg0,arg1)
4356        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4357        if isinstance(arg0,numarray.NumArray):
4358           if isinstance(arg1,Symbol):
4359               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4360           else:
4361               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4362                   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)
4363               arg0_c=arg0.copy()
4364               arg1_c=arg1.copy()
4365               sh0,sh1=arg0.shape,arg1.shape
4366               d0,d1,d01=1,1,1
4367               for i in sh0[axis_offset:]: d0*=i
4368               for i in sh1[axis_offset:]: d1*=i
4369               for i in sh0[:axis_offset]: d01*=i
4370               arg0_c.resize((d01,d0))
4371               arg1_c.resize((d01,d1))
4372               out=numarray.zeros((d0,d1),numarray.Float64)
4373               for i0 in range(d0):
4374                        for i1 in range(d1):
4375                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4376               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4377               return out
4378        elif isinstance(arg0,escript.Data):
4379           if isinstance(arg1,Symbol):
4380               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4381           else:
4382               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4383        else:      
4384           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4385                    
4386    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4387       """
4388       Symbol representing the general tensor product of the transposed of arg0 and arg1
4389       """
4390       def __init__(self,arg0,arg1,axis_offset=0):
4391           """
4392           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4393    
4394           @param arg0: first argument
4395           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4396           @param arg1: second argument
4397           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4398           @raise ValueError: inconsistent dimensions of arguments.
4399           @note: if both arguments have a spatial dimension, they must equal.
4400           """
4401           sh_arg0=pokeShape(arg0)
4402           sh_arg1=pokeShape(arg1)
4403           sh01=sh_arg0[:axis_offset]
4404           sh10=sh_arg1[:axis_offset]
4405           sh0=sh_arg0[axis_offset:]
4406           sh1=sh_arg1[axis_offset:]
4407           if not sh01==sh10:
4408               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)
4409           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4410    
4411       def getMyCode(self,argstrs,format="escript"):
4412          """
4413          returns a program code that can be used to evaluate the symbol.
4414    
4415          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4416          @type argstrs: C{list} of length 2 of C{str}.
4417          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4418          @type format: C{str}
4419          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4420          @rtype: C{str}
4421          @raise NotImplementedError: if the requested format is not available
4422          """
4423          if format=="escript" or format=="str" or format=="text":
4424             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4425          else:
4426             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4427    
4428       def substitute(self,argvals):
4429          """
4430          assigns new values to symbols in the definition of the symbol.
4431          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4432    
4433          @param argvals: new values assigned to symbols
4434          @type argvals: C{dict} with keywords of type L{Symbol}.
4435          @return: result of the substitution process. Operations are executed as much as possible.
4436          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4437          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4438          """
4439          if argvals.has_key(self):
4440             arg=argvals[self]
4441             if self.isAppropriateValue(arg):
4442                return arg
4443             else:
4444                raise TypeError,"%s: new value is not appropriate."%str(self)
4445          else:
4446             args=self.getSubstitutedArguments(argvals)
4447             return generalTransposedTensorProduct(args[0],args[1],args[2])
4448    
4449    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4450        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4451        # calculate the return shape:
4452        shape01=arg0.getShape()[:axis_offset]
4453        shape10=arg1.getShape()[:axis_offset]
4454        shape0=arg0.getShape()[axis_offset:]
4455        shape1=arg1.getShape()[axis_offset:]
4456        if not shape01==shape10:
4457            raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4458    
4459        # whatr function space should be used? (this here is not good!)
4460        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4461        # create return value:
4462        out=escript.Data(0.,tuple(shape0+shape1),fs)
4463        #
4464        s0=[[]]
4465        for k in shape0:
4466              s=[]
4467              for j in s0:
4468                    for i in range(k): s.append(j+[slice(i,i)])
4469              s0=s
4470        s1=[[]]
4471        for k in shape1:
4472              s=[]
4473              for j in s1:
4474                    for i in range(k): s.append(j+[slice(i,i)])
4475              s1=s
4476        s01=[[]]
4477        for k in shape01:
4478              s=[]
4479              for j in s01:
4480                    for i in range(k): s.append(j+[slice(i,i)])
4481              s01=s
4482    
4483        for i0 in s0:
4484           for i1 in s1:
4485             s=escript.Scalar(0.,fs)
4486             for i01 in s01:
4487                s+=arg0.__getitem__(tuple(i01+i0))*arg1.__getitem__(tuple(i01+i1))
4488             out.__setitem__(tuple(i0+i1),s)
4489        return out
4490    
4491    
4492    def matrix_transposed_mult(arg0,arg1):
4493        """
4494        matrix-transposed(matrix) product of the two argument:
4495    
4496        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4497    
4498        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4499    
4500        The last dimensions of arg0 and arg1 must match.
4501    
4502        @param arg0: first argument of rank 2
4503        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4504        @param arg1: second argument of rank 2
4505        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4506        @return: the product of arg0 and the transposed of arg1 at each data point
4507        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4508        @raise ValueError: if the shapes of the arguments are not appropriate
4509        """
4510        sh0=pokeShape(arg0)
4511        sh1=pokeShape(arg1)
4512        if not len(sh0)==2 :
4513            raise ValueError,"first argument must have rank 2"
4514        if not len(sh1)==2 and not len(sh1)==1:
4515            raise ValueError,"second argument must have rank 1 or 2"
4516        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4517    
4518    def tensor_transposed_mult(arg0,arg1):
4519        """
4520        the tensor product of the first and the transpose of the second argument
4521        
4522        for arg0 of rank 2 this is
4523    
4524        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4525    
4526        and for arg0 of rank 4 this is
4527    
4528        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4529                  
4530        or
4531    
4532        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4533    
4534        In the first case the the second dimension of arg0 and arg1 must match and  
4535        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4536    
4537        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4538    
4539        @param arg0: first argument of rank 2 or 4
4540        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4541        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4542        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4543        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4544        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4545        """
4546        sh0=pokeShape(arg0)
4547        sh1=pokeShape(arg1)
4548        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4549           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4550        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4551           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4552        else:
4553            raise ValueError,"first argument must have rank 2 or 4"
4554    
4555    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4556        """
4557        generalized tensor product of transposed of arg0 and arg1:
4558    
4559        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4560    
4561        where
4562    
4563            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4564            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4565            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4566    
4567        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4568        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4569    
4570        @param arg0: first argument
4571        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4572        @param arg1: second argument
4573        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4574        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4575        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4576        """
4577        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4578        arg0,arg1=matchType(arg0,arg1)
4579        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4580        if isinstance(arg0,numarray.NumArray):
4581           if isinstance(arg1,Symbol):
4582               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4583           else:
4584               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4585                   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)
4586               arg0_c=arg0.copy()
4587               arg1_c=arg1.copy()
4588               sh0,sh1=arg0.shape,arg1.shape
4589               d0,d1,d01=1,1,1
4590               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4591               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4592               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4593               arg0_c.resize((d0,d01))
4594               arg1_c.resize((d1,d01))
4595               out=numarray.zeros((d0,d1),numarray.Float64)
4596               for i0 in range(d0):
4597                        for i1 in range(d1):
4598                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4599               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4600               return out
4601        elif isinstance(arg0,escript.Data):
4602           if isinstance(arg1,Symbol):
4603               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4604           else:
4605               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4606        else:      
4607           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4608                    
4609    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4610       """
4611       Symbol representing the general tensor product of arg0 and the transpose of arg1
4612       """
4613       def __init__(self,arg0,arg1,axis_offset=0):
4614           """
4615           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4616    
4617           @param arg0: first argument
4618           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4619           @param arg1: second argument
4620           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4621           @raise ValueError: inconsistent dimensions of arguments.
4622           @note: if both arguments have a spatial dimension, they must equal.
4623           """
4624           sh_arg0=pokeShape(arg0)
4625           sh_arg1=pokeShape(arg1)
4626           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4627           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4628           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4629           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4630           if not sh01==sh10:
4631               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)
4632           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4633    
4634       def getMyCode(self,argstrs,format="escript"):
4635          """
4636          returns a program code that can be used to evaluate the symbol.
4637    
4638          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4639          @type argstrs: C{list} of length 2 of C{str}.
4640          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4641          @type format: C{str}
4642          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4643          @rtype: C{str}
4644          @raise NotImplementedError: if the requested format is not available
4645          """
4646          if format=="escript" or format=="str" or format=="text":
4647             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4648          else:
4649             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4650    
4651       def substitute(self,argvals):
4652          """
4653          assigns new values to symbols in the definition of the symbol.
4654          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4655    
4656          @param argvals: new values assigned to symbols
4657          @type argvals: C{dict} with keywords of type L{Symbol}.
4658          @return: result of the substitution process. Operations are executed as much as possible.
4659          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4660          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4661          """
4662          if argvals.has_key(self):
4663             arg=argvals[self]
4664             if self.isAppropriateValue(arg):
4665                return arg
4666             else:
4667                raise TypeError,"%s: new value is not appropriate."%str(self)
4668          else:
4669             args=self.getSubstitutedArguments(argvals)
4670             return generalTensorTransposedProduct(args[0],args[1],args[2])
4671    
4672    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4673        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4674        # calculate the return shape:
4675        shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4676        shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4677        shape10=arg1.getShape()[arg1.getRank()-axis_offset:]
4678        shape1=arg1.getShape()[:arg1.getRank()-axis_offset]
4679        if not shape01==shape10:
4680            raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4681    
4682        # whatr function space should be used? (this here is not good!)
4683        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4684        # create return value:
4685        out=escript.Data(0.,tuple(shape0+shape1),fs)
4686        #
4687        s0=[[]]
4688        for k in shape0:
4689              s=[]
4690              for j in s0:
4691                    for i in range(k): s.append(j+[slice(i,i)])
4692              s0=s
4693        s1=[[]]
4694        for k in shape1:
4695              s=[]
4696              for j in s1:
4697                    for i in range(k): s.append(j+[slice(i,i)])
4698              s1=s
4699        s01=[[]]
4700        for k in shape01:
4701              s=[]
4702              for j in s01:
4703                    for i in range(k): s.append(j+[slice(i,i)])
4704              s01=s
4705    
4706        for i0 in s0:
4707           for i1 in s1:
4708             s=escript.Scalar(0.,fs)
4709             for i01 in s01:
4710                s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i1+i01))
4711             out.__setitem__(tuple(i0+i1),s)
4712        return out
4713    
4714    
4715  #=========================================================  #=========================================================
4716  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
4717  #=========================================================  #=========================================================
# Line 3776  def grad(arg,where=None): Line 4732  def grad(arg,where=None):
4732                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4733      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4734      @return: gradient of arg.      @return: gradient of arg.
4735      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4736      """      """
4737      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4738         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3804  class Grad_Symbol(DependendSymbol): Line 4760  class Grad_Symbol(DependendSymbol):
4760        d=arg.getDim()        d=arg.getDim()
4761        if d==None:        if d==None:
4762           raise ValueError,"argument must have a spatial dimension"           raise ValueError,"argument must have a spatial dimension"
4763        super(Grad_Symbol,self).__init__(args=[arg,where],shape=tuple(list(arg.getShape()).extend(d)),dim=d)        super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4764    
4765     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4766        """        """
# Line 3816  class Grad_Symbol(DependendSymbol): Line 4772  class Grad_Symbol(DependendSymbol):
4772        @type format: C{str}        @type format: C{str}
4773        @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.
4774        @rtype: C{str}        @rtype: C{str}
4775        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4776        """        """
4777        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4778           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 3869  def integrate(arg,where=None): Line 4825  def integrate(arg,where=None):
4825                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4826      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4827      @return: integral of arg.      @return: integral of arg.
4828      @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}      @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4829      """      """
4830      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4831         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
# Line 3907  class Integrate_Symbol(DependendSymbol): Line 4863  class Integrate_Symbol(DependendSymbol):
4863        @type format: C{str}        @type format: C{str}
4864        @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.
4865        @rtype: C{str}        @rtype: C{str}
4866        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4867        """        """
4868        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4869           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 3959  def interpolate(arg,where): Line 4915  def interpolate(arg,where):
4915      @param where: FunctionSpace to be interpolated to      @param where: FunctionSpace to be interpolated to
4916      @type where: L{escript.FunctionSpace}      @type where: L{escript.FunctionSpace}
4917      @return: interpolated argument      @return: interpolated argument
4918      @rtype:  C{escript.Data} or L{Symbol}      @rtype: C{escript.Data} or L{Symbol}
4919      """      """
4920      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4921         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
# Line 3990  class Interpolate_Symbol(DependendSymbol Line 4946  class Interpolate_Symbol(DependendSymbol
4946        @type format: C{str}        @type format: C{str}
4947        @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.
4948        @rtype: C{str}        @rtype: C{str}
4949        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4950        """        """
4951        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4952           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4043  def div(arg,where=None): Line 4999  def div(arg,where=None):
4999                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
5000      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
5001      @return: divergence of arg.      @return: divergence of arg.
5002      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5003      """      """
5004      if not arg.getShape()==(arg.getDim(),):      if isinstance(arg,Symbol):
5005        raise ValueError,"div: expected shape is (%s,)"%arg.getDim()          dim=arg.getDim()
5006        elif isinstance(arg,escript.Data):
5007            dim=arg.getDomain().getDim()
5008        else:
5009            raise TypeError,"div: argument type not supported"
5010        if not arg.getShape()==(dim,):
5011          raise ValueError,"div: expected shape is (%s,)"%dim
5012      return trace(grad(arg,where))      return trace(grad(arg,where))
5013    
5014  def jump(arg,domain=None):  def jump(arg,domain=None):
# Line 4059  def jump(arg,domain=None): Line 5021  def jump(arg,domain=None):
5021                     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.
5022      @type domain: C{None} or L{escript.Domain}      @type domain: C{None} or L{escript.Domain}
5023      @return: jump of arg      @return: jump of arg
5024      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5025      """      """
5026      if domain==None: domain=arg.getDomain()      if domain==None: domain=arg.getDomain()
5027      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
 #=============================  
 #  
 # wrapper for various functions: if the argument has attribute the function name  
 # as an argument it calls the corresponding methods. Otherwise the corresponding  
 # numarray function is called.  
5028    
5029  # functions involving the underlying Domain:  def L2(arg):
   
 def transpose(arg,axis=None):  
5030      """      """
5031      Returns the transpose of the Data object arg.      returns the L2 norm of arg at where
5032        
5033      @param arg:      @param arg: function which L2 to be calculated.
5034        @type arg: L{escript.Data} or L{Symbol}
5035        @return: L2 norm of arg.
5036        @rtype: L{float} or L{Symbol}
5037        @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5038      """      """
5039      if axis==None:      return sqrt(integrate(inner(arg,arg)))
5040         r=0  #=============================
5041         if hasattr(arg,"getRank"): r=arg.getRank()  #
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
     if isinstance(arg,Symbol):  
        return Transpose_Symbol(arg,axis=r)  
     if isinstance(arg,escript.Data):  
        # hack for transpose  
        r=arg.getRank()  
        if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"  
        s=arg.getShape()  
        out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())  
        for i in range(s[0]):  
           for j in range(s[1]):  
              out[j,i]=arg[i,j]  
        return out  
        # end hack for transpose  
        return arg.transpose(axis)  
     else:  
        return numarray.transpose(arg,axis=axis)  
   
   
5042    
5043  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5044      """      """
5045      resorts the component of arg according to index      resorts the component of arg according to index
5046    
5047      """      """
5048      pass      raise NotImplementedError
5049  #  #
5050  # $Log: util.py,v $  # $Log: util.py,v $
5051  # Revision 1.14.2.16  2005/10/19 06:09:57  gross  # Revision 1.14.2.16  2005/10/19 06:09:57  gross

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

  ViewVC Help
Powered by ViewVC 1.1.26