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

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

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

revision 341 by gross, Mon Dec 12 05:26:10 2005 UTC revision 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: 329 $"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
23    
24    
# 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 maximum(arg0,arg1):  
 # def minimum(arg0,arg1):  
   
 # def transpose(arg,axis=None):  
 # def trace(arg,axis0=0,axis1=1):  
 # def reorderComponents(arg,index):  
   
 # def integrate(arg,where=None):  
 # def interpolate(arg,where):  
 # def div(arg,where=None):  
 # def grad(arg,where=None):  
   
 #  
 # slicing: get  
 #          set  
 #  
 # and derivatives  
   
30  #=========================================================  #=========================================================
31  #   some helpers:  #   some helpers:
32  #=========================================================  #=========================================================
# Line 68  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 96  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 125  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 143  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 163  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 178  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 191  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 249  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 338  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 360  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      try:      if isinstance(arg,numarray.NumArray):
337           return not Lsup(arg)>0.
338        elif isinstance(arg,escript.Data):
339           return False
340        elif isinstance(arg,float):
341           return not Lsup(arg)>0.
342        elif isinstance(arg,int):
343         return not Lsup(arg)>0.         return not Lsup(arg)>0.
344      except TypeError:      elif isinstance(arg,Symbol):
345           return False
346        else:
347         return False         return False
348    
349  def matchType(arg0=0.,arg1=0.):  def matchType(arg0=0.,arg1=0.):
# Line 386  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 412  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 456  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
   
     If shape is not given the shape "largest" shape of args is used.  
438    
439      @param args: a given ob      @param arg0: a given object
440      @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
441      @return: True if the argument is identical to zero.      @param arg1: a given object
442      @rtype: C{list} of C{int}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
443        @return: arg0 and arg1 where copies are returned when the shape has to be changed.
444        @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 491  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 535  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 544  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 568  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 597  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 687  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 695  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 704  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 721  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 732  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 743  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 754  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 765  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 776  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 787  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 798  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 809  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 820  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 875  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 903  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        if arg.rank==0:        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
988           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
989             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))  
990     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
991        return arg._wherePositive()        return arg._wherePositive()
992     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 953  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 989  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        if arg.rank==0:        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1070           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1071             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))  
1072     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1073        return arg._whereNegative()        return arg._whereNegative()
1074     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1039  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 1075  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        if arg.rank==0:        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1152           if arg<0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1153             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))  
1154     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1155        return arg._whereNonNegative()        return arg._whereNonNegative()
1156     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1109  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        if arg.rank==0:        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1182           if arg>0:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1183             return numarray.array(0.)        return out
          else:  
            return numarray.array(1.)  
       else:  
          return numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.  
1184     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1185        return arg._whereNonPositive()        return arg._whereNonPositive()
1186     elif isinstance(arg,float):     elif isinstance(arg,float):
# Line 1145  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        if arg.rank==0:        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1214           if abs(arg)<=tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1215             return numarray.array(1.)        return out
          else:  
            return numarray.array(0.)  
       else:  
          return numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
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 1198  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 1232  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        if arg.rank==0:        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1294          if abs(arg)>tol:        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1295             return numarray.array(1.)        return out
         else:  
            return numarray.array(0.)  
       else:  
          return numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.  
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 1269  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 1307  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 1359  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 1397  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 1449  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 1487  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 1539  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 1577  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 1629  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 1667  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 1719  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 1757  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 1809  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 1847  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 1899  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 1937  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 1989  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 2027  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 2079  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 2117  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 2169  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 2207  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 2259  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 2297  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 2349  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 2387  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 2439  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 2477  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 2529  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 2567  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 2619  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 2667  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 2719  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 2760  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 2796  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 2837  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 2873  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    
2927    def trace(arg,axis_offset=0):
2928       """
2929       returns the trace of arg which the sum of arg[k,k] over k.
2930    
2931       @param arg: argument
2932       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2933       @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
2934                      axis_offset and axis_offset+1 must be equal.
2935       @type axis_offset: C{int}
2936       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2937       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2938       """
2939       if isinstance(arg,numarray.NumArray):
2940          sh=arg.shape
2941          if len(sh)<2:
2942            raise ValueError,"rank of argument must be greater than 1"
2943          if axis_offset<0 or axis_offset>len(sh)-2:
2944            raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
2945          s1=1
2946          for i in range(axis_offset): s1*=sh[i]
2947          s2=1
2948          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2949          if not sh[axis_offset] == sh[axis_offset+1]:
2950            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))
2952          out=numarray.zeros([s1,s2],numarray.Float64)
2953          for i1 in range(s1):
2954            for i2 in range(s2):
2955                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:])
2957          return out
2958       elif isinstance(arg,escript.Data):
2959          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):
2968          raise TypeError,"illegal argument type float."
2969       elif isinstance(arg,int):
2970          raise TypeError,"illegal argument type int."
2971       elif isinstance(arg,Symbol):
2972          return Trace_Symbol(arg,axis_offset)
2973       else:
2974          raise TypeError,"Unknown argument type."
2975    
2976    class Trace_Symbol(DependendSymbol):
2977       """
2978       L{Symbol} representing the result of the trace function
2979       """
2980       def __init__(self,arg,axis_offset=0):
2981          """
2982          initialization of trace L{Symbol} with argument arg
2983          @param arg: argument of function
2984          @type arg: L{Symbol}.
2985          @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
2986                      axis_offset and axis_offset+1 must be equal.
2987          @type axis_offset: C{int}
2988          """
2989          if arg.getRank()<2:
2990            raise ValueError,"rank of argument must be greater than 1"
2991          if axis_offset<0 or axis_offset>arg.getRank()-2:
2992            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
2993          s=list(arg.getShape())        
2994          if not s[axis_offset] == s[axis_offset+1]:
2995            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())
2997    
2998       def getMyCode(self,argstrs,format="escript"):
2999          """
3000          returns a program code that can be used to evaluate the symbol.
3001    
3002          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3003          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3004          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3005          @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.
3007          @rtype: C{str}
3008          @raise NotImplementedError: if the requested format is not available
3009          """
3010          if format=="escript" or format=="str"  or format=="text":
3011             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3012          else:
3013             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
3014    
3015       def substitute(self,argvals):
3016          """
3017          assigns new values to symbols in the definition of the symbol.
3018          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3019    
3020          @param argvals: new values assigned to symbols
3021          @type argvals: C{dict} with keywords of type L{Symbol}.
3022          @return: result of the substitution process. Operations are executed as much as possible.
3023          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3024          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3025          """
3026          if argvals.has_key(self):
3027             arg=argvals[self]
3028             if self.isAppropriateValue(arg):
3029                return arg
3030             else:
3031                raise TypeError,"%s: new value is not appropriate."%str(self)
3032          else:
3033             arg=self.getSubstitutedArguments(argvals)
3034             return trace(arg[0],axis_offset=arg[1])
3035    
3036       def diff(self,arg):
3037          """
3038          differential of this object
3039    
3040          @param arg: the derivative is calculated with respect to arg
3041          @type arg: L{escript.Symbol}
3042          @return: derivative with respect to C{arg}
3043          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3044          """
3045          if arg==self:
3046             return identity(self.getShape())
3047          else:
3048             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 2936  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 2995  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 3035  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 3095  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 3140  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 3201  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 3244  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 3304  def maximum(*args): Line 3912  def maximum(*args):
3912         if out==None:         if out==None:
3913            out=a            out=a
3914         else:         else:
3915            m=whereNegative(out-a)            diff=add(a,-out)
3916            out=m*a+(1.-m)*out            out=add(out,mult(wherePositive(diff),diff))
3917      return out      return out
3918        
3919  def minimum(*arg):  def minimum(*args):
3920      """      """
3921      the minimum over arguments args      the minimum over arguments args
3922    
# Line 3322  def minimum(*arg): Line 3930  def minimum(*arg):
3930         if out==None:         if out==None:
3931            out=a            out=a
3932         else:         else:
3933            m=whereNegative(out-a)            diff=add(a,-out)
3934            out=m*out+(1.-m)*a            out=add(out,mult(whereNegative(diff),diff))
3935      return out      return out
3936    
3937    def clip(arg,minval=0.,maxval=1.):
3938        """
3939        cuts the values of arg between minval and maxval
3940    
3941        @param arg: argument
3942        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
3943        @param minval: lower range
3944        @type minval: C{float}
3945        @param maxval: upper range
3946        @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
3948                 less then maxval are unchanged.
3949        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
3950        @raise ValueError: if minval>maxval
3951        """
3952        if minval>maxval:
3953           raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
3954        return minimum(maximum(minval,arg),maxval)
3955    
3956        
3957  def inner(arg0,arg1):  def inner(arg0,arg1):
3958      """      """
# Line 3340  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 3348  def inner(arg0,arg1): Line 3976  def inner(arg0,arg1):
3976      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
3977      if not sh0==sh1:      if not sh0==sh1:
3978          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
3979      return generalTensorProduct(arg0,arg1,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 3376  def matrixmult(arg0,arg1): Line 4030  def matrixmult(arg0,arg1):
4030          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4031      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
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,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,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 3415  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 3436  def tensormult(arg0,arg1): Line 4076  def tensormult(arg0,arg1):
4076      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
4077      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
4078      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4079         return generalTensorProduct(arg0,arg1,offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
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,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,offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4086      """      """
4087      generalized tensor product      generalized tensor product
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-offset]      where
           r runs trough arg0.Shape[:offset]  
           t runs through arg1.Shape[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 3467  def generalTensorProduct(arg0,arg1,offse Line 4106  def generalTensorProduct(arg0,arg1,offse
4106      # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols      # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4107      if isinstance(arg0,numarray.NumArray):      if isinstance(arg0,numarray.NumArray):
4108         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4109             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4110         else:         else:
4111             if not arg0.shape[arg0.rank-offset:]==arg1.shape[: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."%(offset,offset)                 raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
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
4116             d0,d1,d01=1,1,1             d0,d1,d01=1,1,1
4117             for i in sh0[:arg0.rank-offset]: d0*=i             for i in sh0[:arg0.rank-axis_offset]: d0*=i
4118             for i in sh1[offset:]: d1*=i             for i in sh1[axis_offset:]: d1*=i
4119             for i in sh1[: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])
4126             out.resize(sh0[:arg0.rank-offset]+sh1[offset:])             out.resize(sh0[:arg0.rank-axis_offset]+sh1[axis_offset:])
4127             return out             return out
4128      elif isinstance(arg0,escript.Data):      elif isinstance(arg0,escript.Data):
4129         if isinstance(arg1,Symbol):         if isinstance(arg1,Symbol):
4130             return GeneralTensorProduct_Symbol(arg0,arg1,offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4131         else:         else:
4132             return escript_generalTensorProduct(arg0,arg1,offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,offset)             return escript_generalTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4133      else:            else:      
4134         return GeneralTensorProduct_Symbol(arg0,arg1,offset)         return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
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,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)
4152         sh_arg1=pokeShape(arg1)         sh_arg1=pokeShape(arg1)
4153         sh0=sh_arg0[:len(sh_arg0)-offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4154         sh01=sh_arg0[len(sh_arg0)-offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4155         sh10=sh_arg1[:offset]         sh10=sh_arg1[:axis_offset]
4156         sh1=sh_arg1[offset:]         sh1=sh_arg1[axis_offset:]
4157         if not sh01==sh10:         if not sh01==sh10:
4158             raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)             raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4159         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,offset])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4160    
4161     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
4162        """        """
# Line 3529  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,offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])           return "generalTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4175        else:        else:
4176           raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)           raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4177    
# Line 3557  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_generalTensorProduct(arg0,arg1,offset): # this should be escript._generalTensorProduct  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
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:
4216      shape0=arg0.getShape()[:arg0.getRank()-offset]      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4217      shape01=arg0.getShape()[arg0.getRank()-offset:]      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4218      shape10=arg1.getShape()[:offset]      shape10=arg1.getShape()[:axis_offset]
4219      shape1=arg1.getShape()[offset:]      shape1=arg1.getShape()[axis_offset:]
4220      if not shape01==shape10:      if not shape01==shape10:
4221          raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(offset,offset)          raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4222    
4223        # whatr function space should be used? (this here is not good!)
4224        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4225      # create return value:      # create return value:
4226      out=escript.Data(0.,tuple(shape0+shape1),arg0.getFunctionSpace())      out=escript.Data(0.,tuple(shape0+shape1),fs)
4227      #      #
4228      s0=[[]]      s0=[[]]
4229      for k in shape0:      for k in shape0:
# Line 3591  def escript_generalTensorProduct(arg0,ar Line 4246  def escript_generalTensorProduct(arg0,ar
4246    
4247      for i0 in s0:      for i0 in s0:
4248         for i1 in s1:         for i1 in s1:
4249           s=escript.Scalar(0.,arg0.getFunctionSpace())           s=escript.Scalar(0.,fs)
4250           for i01 in s01:           for i01 in s01:
4251              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))
4252           out.__setitem__(tuple(i0+i1),s)           out.__setitem__(tuple(i0+i1),s)
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  #   some little helpers  #  functions dealing with spatial dependency
4717  #=========================================================  #=========================================================
4718  def grad(arg,where=None):  def grad(arg,where=None):
4719      """      """
4720      Returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
4721    
4722        If C{g} is the returned object, then
4723    
4724          - if C{arg} is rank 0 C{g[s]} is the derivative of C{arg} with respect to the C{s}-th spatial dimension.
4725          - if C{arg} is rank 1 C{g[i,s]} is the derivative of C{arg[i]} with respect to the C{s}-th spatial dimension.
4726          - if C{arg} is rank 2 C{g[i,j,s]} is the derivative of C{arg[i,j]} with respect to the C{s}-th spatial dimension.
4727          - if C{arg} is rank 3 C{g[i,j,k,s]} is the derivative of C{arg[i,j,k]} with respect to the C{s}-th spatial dimension.
4728    
4729      @param arg:   Data object representing the function which gradient      @param arg: function which gradient to be calculated. Its rank has to be less than 3.
4730                    to be calculated.      @type arg: L{escript.Data} or L{Symbol}
4731      @param where: FunctionSpace in which the gradient will be calculated.      @param where: FunctionSpace in which the gradient will be calculated.
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}
4734        @return: gradient of arg.
4735        @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 3617  def grad(arg,where=None): Line 4742  def grad(arg,where=None):
4742         else:         else:
4743            return arg._grad(where)            return arg._grad(where)
4744      else:      else:
4745        raise TypeError,"grad: Unknown argument type."         raise TypeError,"grad: Unknown argument type."
4746    
4747    class Grad_Symbol(DependendSymbol):
4748       """
4749       L{Symbol} representing the result of the gradient operator
4750       """
4751       def __init__(self,arg,where=None):
4752          """
4753          initialization of gradient L{Symbol} with argument arg
4754          @param arg: argument of function
4755          @type arg: L{Symbol}.
4756          @param where: FunctionSpace in which the gradient will be calculated.
4757                      If not present or C{None} an appropriate default is used.
4758          @type where: C{None} or L{escript.FunctionSpace}
4759          """
4760          d=arg.getDim()
4761          if d==None:
4762             raise ValueError,"argument must have a spatial dimension"
4763          super(Grad_Symbol,self).__init__(args=[arg,where],shape=arg.getShape()+(d,),dim=d)
4764    
4765       def getMyCode(self,argstrs,format="escript"):
4766          """
4767          returns a program code that can be used to evaluate the symbol.
4768    
4769          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4770          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4771          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4772          @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.
4774          @rtype: C{str}
4775          @raise NotImplementedError: if the requested format is not available
4776          """
4777          if format=="escript" or format=="str"  or format=="text":
4778             return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
4779          else:
4780             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4781    
4782       def substitute(self,argvals):
4783          """
4784          assigns new values to symbols in the definition of the symbol.
4785          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4786    
4787          @param argvals: new values assigned to symbols
4788          @type argvals: C{dict} with keywords of type L{Symbol}.
4789          @return: result of the substitution process. Operations are executed as much as possible.
4790          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4791          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4792          """
4793          if argvals.has_key(self):
4794             arg=argvals[self]
4795             if self.isAppropriateValue(arg):
4796                return arg
4797             else:
4798                raise TypeError,"%s: new value is not appropriate."%str(self)
4799          else:
4800             arg=self.getSubstitutedArguments(argvals)
4801             return grad(arg[0],where=arg[1])
4802    
4803       def diff(self,arg):
4804          """
4805          differential of this object
4806    
4807          @param arg: the derivative is calculated with respect to arg
4808          @type arg: L{escript.Symbol}
4809          @return: derivative with respect to C{arg}
4810          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4811          """
4812          if arg==self:
4813             return identity(self.getShape())
4814          else:
4815             return grad(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4816    
4817  def integrate(arg,where=None):  def integrate(arg,where=None):
4818      """      """
4819      Return the integral if the function represented by Data object arg over      Return the integral of the function C{arg} over its domain. If C{where} is present C{arg} is interpolated to C{where}
4820      its domain.      before integration.
4821    
4822      @param arg:   Data object representing the function which is integrated.      @param arg:   the function which is integrated.
4823        @type arg: L{escript.Data} or L{Symbol}
4824      @param where: FunctionSpace in which the integral is calculated.      @param where: FunctionSpace in which the integral is calculated.
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}
4827        @return: integral of arg.
4828        @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4829      """      """
4830      if isinstance(arg,numarray.NumArray):      if isinstance(arg,Symbol):
         if checkForZero(arg):  
            return arg  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,float):  
         if checkForZero(arg):  
            return arg  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,int):  
         if checkForZero(arg):  
            return float(arg)  
         else:  
            raise TypeError,"integrate: cannot intergrate argument"  
     elif isinstance(arg,Symbol):  
4831         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
4832      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
4833         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
# Line 3654  def integrate(arg,where=None): Line 4838  def integrate(arg,where=None):
4838      else:      else:
4839        raise TypeError,"integrate: Unknown argument type."        raise TypeError,"integrate: Unknown argument type."
4840    
4841    class Integrate_Symbol(DependendSymbol):
4842       """
4843       L{Symbol} representing the result of the spatial integration operator
4844       """
4845       def __init__(self,arg,where=None):
4846          """
4847          initialization of integration L{Symbol} with argument arg
4848          @param arg: argument of the integration
4849          @type arg: L{Symbol}.
4850          @param where: FunctionSpace in which the integration will be calculated.
4851                      If not present or C{None} an appropriate default is used.
4852          @type where: C{None} or L{escript.FunctionSpace}
4853          """
4854          super(Integrate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4855    
4856       def getMyCode(self,argstrs,format="escript"):
4857          """
4858          returns a program code that can be used to evaluate the symbol.
4859    
4860          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4861          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4862          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4863          @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.
4865          @rtype: C{str}
4866          @raise NotImplementedError: if the requested format is not available
4867          """
4868          if format=="escript" or format=="str"  or format=="text":
4869             return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
4870          else:
4871             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4872    
4873       def substitute(self,argvals):
4874          """
4875          assigns new values to symbols in the definition of the symbol.
4876          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4877    
4878          @param argvals: new values assigned to symbols
4879          @type argvals: C{dict} with keywords of type L{Symbol}.
4880          @return: result of the substitution process. Operations are executed as much as possible.
4881          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4882          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4883          """
4884          if argvals.has_key(self):
4885             arg=argvals[self]
4886             if self.isAppropriateValue(arg):
4887                return arg
4888             else:
4889                raise TypeError,"%s: new value is not appropriate."%str(self)
4890          else:
4891             arg=self.getSubstitutedArguments(argvals)
4892             return integrate(arg[0],where=arg[1])
4893    
4894       def diff(self,arg):
4895          """
4896          differential of this object
4897    
4898          @param arg: the derivative is calculated with respect to arg
4899          @type arg: L{escript.Symbol}
4900          @return: derivative with respect to C{arg}
4901          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
4902          """
4903          if arg==self:
4904             return identity(self.getShape())
4905          else:
4906             return integrate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4907    
4908    
4909  def interpolate(arg,where):  def interpolate(arg,where):
4910      """      """
4911      Interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where.
4912    
4913      @param arg:    interpolant      @param arg: interpolant
4914      @param where:  FunctionSpace to interpolate to      @type arg: L{escript.Data} or L{Symbol}
4915        @param where: FunctionSpace to be interpolated to
4916        @type where: L{escript.FunctionSpace}
4917        @return: interpolated argument
4918        @rtype: C{escript.Data} or L{Symbol}
4919      """      """
4920      if testForZero(arg):      if isinstance(arg,Symbol):
4921        return 0         return Interpolate_Symbol(arg,where)
     elif isinstance(arg,Symbol):  
        return Interpolated_Symbol(arg,where)  
4922      else:      else:
4923         return escript.Data(arg,where)         return escript.Data(arg,where)
4924    
4925  def div(arg,where=None):  class Interpolate_Symbol(DependendSymbol):
4926      """     """
4927      Returns the divergence of arg at where.     L{Symbol} representing the result of the interpolation operator
4928       """
4929       def __init__(self,arg,where):
4930          """
4931          initialization of interpolation L{Symbol} with argument arg
4932          @param arg: argument of the interpolation
4933          @type arg: L{Symbol}.
4934          @param where: FunctionSpace into which the argument is interpolated.
4935          @type where: L{escript.FunctionSpace}
4936          """
4937          super(Interpolate_Symbol,self).__init__(args=[arg,where],shape=arg.getShape(),dim=arg.getDim())
4938    
4939      @param arg:   Data object representing the function which gradient to     def getMyCode(self,argstrs,format="escript"):
4940                    be calculated.        """
4941      @param where: FunctionSpace in which the gradient will be calculated.        returns a program code that can be used to evaluate the symbol.
                   If not present or C{None} an appropriate default is used.  
     """  
     g=grad(arg,where)  
     return trace(g,axis0=g.getRank()-2,axis1=g.getRank()-1)  
4942    
4943  def jump(arg):        @param argstrs: gives for each argument a string representing the argument for the evaluation.
4944      """        @type argstrs: C{str} or a C{list} of length 1 of C{str}.
4945      Returns the jump of arg across a continuity.        @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
4946          @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.
4948          @rtype: C{str}
4949          @raise NotImplementedError: if the requested format is not available
4950          """
4951          if format=="escript" or format=="str"  or format=="text":
4952             return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
4953          else:
4954             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
4955    
4956      @param arg:   Data object representing the function which gradient     def substitute(self,argvals):
4957                    to be calculated.        """
4958      """        assigns new values to symbols in the definition of the symbol.
4959      d=arg.getDomain()        The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
     return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())  
4960    
4961  #=============================        @param argvals: new values assigned to symbols
4962  #        @type argvals: C{dict} with keywords of type L{Symbol}.
4963  # wrapper for various functions: if the argument has attribute the function name        @return: result of the substitution process. Operations are executed as much as possible.
4964  # as an argument it calls the corresponding methods. Otherwise the corresponding        @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4965  # numarray function is called.        @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4966          """
4967          if argvals.has_key(self):
4968             arg=argvals[self]
4969             if self.isAppropriateValue(arg):
4970                return arg
4971             else:
4972                raise TypeError,"%s: new value is not appropriate."%str(self)
4973          else:
4974             arg=self.getSubstitutedArguments(argvals)
4975             return interpolate(arg[0],where=arg[1])
4976    
4977       def diff(self,arg):
4978          """
4979          differential of this object
4980    
4981          @param arg: the derivative is calculated with respect to arg
4982          @type arg: L{escript.Symbol}
4983          @return: derivative with respect to C{arg}
4984          @rtype: L{Symbol} but other types such as L{escript.Data}, L{numarray.NumArray}  are possible.
4985          """
4986          if arg==self:
4987             return identity(self.getShape())
4988          else:
4989             return interpolate(self.getDifferentiatedArguments(arg)[0],where=self.getArgument()[1])
4990    
 # functions involving the underlying Domain:  
4991    
4992  def transpose(arg,axis=None):  def div(arg,where=None):
4993      """      """
4994      Returns the transpose of the Data object arg.      returns the divergence of arg at where.
4995    
4996      @param arg:      @param arg: function which divergence to be calculated. Its shape has to be (d,) where d is the spatial dimension.
4997        @type arg: L{escript.Data} or L{Symbol}
4998        @param where: FunctionSpace in which the divergence will be calculated.
4999                      If not present or C{None} an appropriate default is used.
5000        @type where: C{None} or L{escript.FunctionSpace}
5001        @return: divergence of arg.
5002        @rtype: L{escript.Data} or L{Symbol}
5003      """      """
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
5004      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5005         return Transpose_Symbol(arg,axis=r)          dim=arg.getDim()
5006      if isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
5007         # hack for transpose          dim=arg.getDomain().getDim()
        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)  
5008      else:      else:
5009         return numarray.transpose(arg,axis=axis)          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))
5013    
5014  def trace(arg,axis0=0,axis1=1):  def jump(arg,domain=None):
5015      """      """
5016      Return      returns the jump of arg across the continuity of the domain
5017    
5018      @param arg:      @param arg: argument
5019        @type arg: L{escript.Data} or L{Symbol}
5020        @param domain: the domain where the discontinuity is located. If domain is not present or equal to C{None}
5021                       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}
5023        @return: jump of arg
5024        @rtype: L{escript.Data} or L{Symbol}
5025      """      """
5026      if isinstance(arg,Symbol):      if domain==None: domain=arg.getDomain()
5027         s=list(arg.getShape())              return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
        s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])  
        return Trace_Symbol(arg,axis0=axis0,axis1=axis1)  
     elif isinstance(arg,escript.Data):  
        # hack for trace  
        s=arg.getShape()  
        if s[axis0]!=s[axis1]:  
            raise ValueError,"illegal axis in trace"  
        out=escript.Scalar(0.,arg.getFunctionSpace())  
        for i in range(s[axis0]):  
           out+=arg[i,i]  
        return out  
        # end hack for trace  
     else:  
        return numarray.trace(arg,axis0=axis0,axis1=axis1)  
5028    
5029    def L2(arg):
5030        """
5031        returns the L2 norm of arg at where
5032        
5033        @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        return sqrt(integrate(inner(arg,arg)))
5040    #=============================
5041    #
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.341  
changed lines
  Added in v.795

  ViewVC Help
Powered by ViewVC 1.1.26