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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26