/[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 443 by gross, Fri Jan 20 06:22:38 2006 UTC revision 1247 by ksteube, Tue Aug 14 01:29:20 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
 #  
 #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  
 #  
 #   This software is the property of ACcESS.  No part of this code  
 #   may be copied in any form or by any means without the expressed written  
 #   consent of ACcESS.  Copying, use or modification of this software  
 #   by any unauthorised person is illegal unless that  
 #   person has a software license agreement with ACcESS.  
 #  
2    
3  """  """
4  Utility functions for escript  Utility functions for escript
5    
 @remark:  This module is under construction and is still tested!!!  
   
6  @var __author__: name of author  @var __author__: name of author
7  @var __licence__: licence agreement  @var __copyright__: copyrights
8    @var __license__: licence agreement
9  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
10  @var __version__: version  @var __version__: version
11  @var __date__: date of the version  @var __date__: date of the version
12  """  """
13                                                                                                                                                                                                                                                                                                                                                                                                            
14  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
15  __licence__="contact: esys@access.uq.edu.au"  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
16                        http://www.access.edu.au
17                    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19                 http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.iservo.edu.au/esys/escript"  __url__="http://www.iservo.edu.au/esys/escript"
21  __version__="$Revision$"  __version__="$Revision$"
22  __date__="$Date$"  __date__="$Date$"
# Line 30  __date__="$Date$" Line 24  __date__="$Date$"
24    
25  import math  import math
26  import numarray  import numarray
 import numarray.linear_algebra  
27  import escript  import escript
28  import os  import os
29    from esys.escript import C_GeneralTensorProduct
30    from esys.escript import getVersion
31    
32    #=========================================================
33    #   some helpers:
34    #=========================================================
35    def getTagNames(domain):
36        """
37        returns a list of the tag names used by the domain
38    
39  # missing tests:      
40        @param domain: a domain object
41        @type domain: L{escript.Domain}
42        @return: a list of the tag name used by the domain.
43        @rtype: C{list} of C{str}
44        """
45        return [n.strip() for n in domain.showTagNames().split(",") ]
46    
47  # def pokeShape(arg):  def insertTagNames(domain,**kwargs):
48  # def pokeDim(arg):      """
49  # def commonShape(arg0,arg1):      inserts tag names into the domain
 # def commonDim(*args):  
 # def testForZero(arg):  
 # def matchType(arg0=0.,arg1=0.):  
 # def matchShape(arg0,arg1):  
50    
51  # def transpose(arg,axis=None):      @param domain: a domain object
52  # def reorderComponents(arg,index):      @type domain: C{escript.Domain}
53        @keyword <tag name>: tag key assigned to <tag name>
54        @type <tag name>: C{int}
55        """
56        for  k in kwargs:
57             domain.setTagMap(k,kwargs[k])
58    
59  #  def insertTaggedValues(target,**kwargs):
60  # slicing: get      """
61  #          set      inserts tagged values into the tagged using tag names
 #  
 # and derivatives  
62    
63  #=========================================================      @param target: data to be filled by tagged values
64  #   some helpers:      @type target: L{escript.Data}
65  #=========================================================      @keyword <tag name>: value to be used for <tag name>
66        @type <tag name>: C{float} or {numarray.NumArray}
67        @return: C{target}
68        @rtype: L{escript.Data}
69        """
70        for k in kwargs:
71            target.setTaggedValue(k,kwargs[k])
72        return target
73    
74        
75  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
76      """      """
77      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.
78    
79      Example:      Example::
80    
81         tmp=Scalar(..)         tmp=Scalar(..)
82         v=Vector(..)         v=Vector(..)
# Line 88  def saveDX(filename,domain=None,**data): Line 104  def saveDX(filename,domain=None,**data):
104      """      """
105      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.
106    
107      Example:      Example::
108    
109         tmp=Scalar(..)         tmp=Scalar(..)
110         v=Vector(..)         v=Vector(..)
# Line 119  def kronecker(d=3): Line 135  def kronecker(d=3):
135     @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
136     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
137     @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
138     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2.
139     """     """
140     return identityTensor(d)     return identityTensor(d)
141    
# Line 134  def identity(shape=()): Line 150  def identity(shape=()):
150     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
151     """     """
152     if len(shape)>0:     if len(shape)>0:
153        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
154        if len(shape)==1:        if len(shape)==1:
155            for i0 in range(shape[0]):            for i0 in range(shape[0]):
156               out[i0,i0]=1.               out[i0,i0]=1.
# Line 155  def identityTensor(d=3): Line 171  def identityTensor(d=3):
171     @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
172     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
173     @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
174     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 2
175     """     """
176     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
177         return escript.Data(identity((d.getDim(),)),d)         return escript.Data(identity((d.getDim(),)),d)
# Line 171  def identityTensor4(d=3): Line 187  def identityTensor4(d=3):
187     @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
188     @type d: C{int} or any object with a C{getDim} method     @type d: C{int} or any object with a C{getDim} method
189     @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
190     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 4.     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 4.
191     """     """
192     if isinstance(d,escript.FunctionSpace):     if isinstance(d,escript.FunctionSpace):
193         return escript.Data(identity((d.getDim(),d.getDim())),d)         return escript.Data(identity((d.getDim(),d.getDim())),d)
# Line 189  def unitVector(i=0,d=3): Line 205  def unitVector(i=0,d=3):
205     @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
206     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}     @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
207     @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
208     @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 1     @rtype: L{numarray.NumArray} or L{escript.Data} of rank 1
209     """     """
210     return kronecker(d)[i]     return kronecker(d)[i]
211    
# Line 245  def inf(arg): Line 261  def inf(arg):
261    
262      @param arg: argument      @param arg: argument
263      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.      @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
264      @return : minimum value of arg over all components and all data points      @return: minimum value of arg over all components and all data points
265      @rtype: C{float}      @rtype: C{float}
266      @raise TypeError: if type of arg cannot be processed      @raise TypeError: if type of arg cannot be processed
267      """      """
# Line 264  def inf(arg): Line 280  def inf(arg):
280  #=========================================================================  #=========================================================================
281  #   some little helpers  #   some little helpers
282  #=========================================================================  #=========================================================================
283  def pokeShape(arg):  def getRank(arg):
284        """
285        identifies the rank of its argument
286    
287        @param arg: a given object
288        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
289        @return: the rank of the argument
290        @rtype: C{int}
291        @raise TypeError: if type of arg cannot be processed
292        """
293    
294        if isinstance(arg,numarray.NumArray):
295            return arg.rank
296        elif isinstance(arg,escript.Data):
297            return arg.getRank()
298        elif isinstance(arg,float):
299            return 0
300        elif isinstance(arg,int):
301            return 0
302        elif isinstance(arg,Symbol):
303            return arg.getRank()
304        else:
305          raise TypeError,"getShape: cannot identify shape"
306    def getShape(arg):
307      """      """
308      identifies the shape of its argument      identifies the shape of its argument
309    
# Line 286  def pokeShape(arg): Line 325  def pokeShape(arg):
325      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
326          return arg.getShape()          return arg.getShape()
327      else:      else:
328        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
329    
330  def pokeDim(arg):  def pokeDim(arg):
331      """      """
# Line 309  def commonShape(arg0,arg1): Line 348  def commonShape(arg0,arg1):
348      """      """
349      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.
350    
351      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
352      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
353      @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.
354      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
355      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
356      """      """
357      sh0=pokeShape(arg0)      sh0=getShape(arg0)
358      sh1=pokeShape(arg1)      sh1=getShape(arg1)
359      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
360         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
361               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 334  def commonDim(*args): Line 373  def commonDim(*args):
373      """      """
374      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.
375    
376      @param *args: given objects      @param args: given objects
377      @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
378               a spatial dimension C{None} is returned.               a spatial dimension C{None} is returned.
379      @rtype: C{int} or C{None}      @rtype: C{int} or C{None}
# Line 356  def testForZero(arg): Line 395  def testForZero(arg):
395    
396      @param arg: a given object      @param arg: a given object
397      @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}
398      @return : True if the argument is identical to zero.      @return: True if the argument is identical to zero.
399      @rtype : C{bool}      @rtype: C{bool}
400      """      """
401      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
402         return not Lsup(arg)>0.         return not Lsup(arg)>0.
# Line 390  def matchType(arg0=0.,arg1=0.): Line 429  def matchType(arg0=0.,arg1=0.):
429         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
430            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
431         elif isinstance(arg1,float):         elif isinstance(arg1,float):
432            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
433         elif isinstance(arg1,int):         elif isinstance(arg1,int):
434            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
435         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
436            pass            pass
437         else:         else:
# Line 416  def matchType(arg0=0.,arg1=0.): Line 455  def matchType(arg0=0.,arg1=0.):
455         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
456            pass            pass
457         elif isinstance(arg1,float):         elif isinstance(arg1,float):
458            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
459         elif isinstance(arg1,int):         elif isinstance(arg1,int):
460            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
461         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
462            pass            pass
463         else:         else:
464            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
465      elif isinstance(arg0,float):      elif isinstance(arg0,float):
466         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
467            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
468         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
469            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
470         elif isinstance(arg1,float):         elif isinstance(arg1,float):
471            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
472            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
473         elif isinstance(arg1,int):         elif isinstance(arg1,int):
474            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
475            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
476         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
477            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
478         else:         else:
479            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
480      elif isinstance(arg0,int):      elif isinstance(arg0,int):
481         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
482            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
483         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
484            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
485         elif isinstance(arg1,float):         elif isinstance(arg1,float):
486            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
487            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
488         elif isinstance(arg1,int):         elif isinstance(arg1,int):
489            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
490            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
491         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
492            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
493         else:         else:
494            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
495      else:      else:
# Line 460  def matchType(arg0=0.,arg1=0.): Line 499  def matchType(arg0=0.,arg1=0.):
499    
500  def matchShape(arg0,arg1):  def matchShape(arg0,arg1):
501      """      """
502            return representations of arg0 amd arg1 which ahve the same shape
503    
504      If shape is not given the shape "largest" shape of args is used.      @param arg0: a given object
505        @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
506      @param args: a given ob      @param arg1: a given object
507      @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}
508      @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.
509      @rtype: C{list} of C{int}      @rtype: C{tuple}
510      """      """
511      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
512      sh0=pokeShape(arg0)      sh0=getShape(arg0)
513      sh1=pokeShape(arg1)      sh1=getShape(arg1)
514      if len(sh0)<len(sh):      if len(sh0)<len(sh):
515         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
516      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
517         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
518      else:      else:
519         return arg0,arg1         return arg0,arg1
520  #=========================================================  #=========================================================
# Line 495  class Symbol(object): Line 534  class Symbol(object):
534         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
535         symbols or any other object.         symbols or any other object.
536    
537         @param arg: the arguments of the symbol.         @param args: the arguments of the symbol.
538         @type arg: C{list}         @type args: C{list}
539         @param shape: the shape         @param shape: the shape
540         @type shape: C{tuple} of C{int}         @type shape: C{tuple} of C{int}
541         @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 539  class Symbol(object): Line 578  class Symbol(object):
578         """         """
579         the shape of the symbol.         the shape of the symbol.
580    
581         @return : the shape of the symbol.         @return: the shape of the symbol.
582         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
583         """         """
584         return self.__shape         return self.__shape
# Line 548  class Symbol(object): Line 587  class Symbol(object):
587         """         """
588         the spatial dimension         the spatial dimension
589    
590         @return : the spatial dimension         @return: the spatial dimension
591         @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.
592         """         """
593         return self.__dim         return self.__dim
# Line 572  class Symbol(object): Line 611  class Symbol(object):
611         """         """
612         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.
613    
614         @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].
615         @type argvals: C{dict} with keywords of type L{Symbol}.         @type argvals: C{dict} with keywords of type L{Symbol}.
616         @rtype: C{list} of objects         @rtype: C{list} of objects
617         @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.
618         """         """
619         out=[]         out=[]
620         for a in self.getArgument():         for a in self.getArgument():
# Line 599  class Symbol(object): Line 638  class Symbol(object):
638            if isinstance(a,Symbol):            if isinstance(a,Symbol):
639               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
640            else:            else:
641                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
642                if len(s)>0:                if len(s)>0:
643                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
644                else:                else:
645                   out.append(a)                   out.append(a)
646         return out         return out
# Line 691  class Symbol(object): Line 730  class Symbol(object):
730         else:         else:
731            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
732            if len(s)>0:            if len(s)>0:
733               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
734            else:            else:
735               return 0.               return 0.
736    
# Line 699  class Symbol(object): Line 738  class Symbol(object):
738         """         """
739         returns -self.         returns -self.
740    
741         @return:  a S{Symbol} representing the negative of the object         @return:  a L{Symbol} representing the negative of the object
742         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
743         """         """
744         return self*(-1.)         return self*(-1.)
# Line 708  class Symbol(object): Line 747  class Symbol(object):
747         """         """
748         returns +self.         returns +self.
749    
750         @return:  a S{Symbol} representing the positive of the object         @return:  a L{Symbol} representing the positive of the object
751         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
752         """         """
753         return self*(1.)         return self*(1.)
754    
755     def __abs__(self):     def __abs__(self):
756         """         """
757         returns a S{Symbol} representing the absolute value of the object.         returns a L{Symbol} representing the absolute value of the object.
758         """         """
759         return Abs_Symbol(self)         return Abs_Symbol(self)
760    
# Line 725  class Symbol(object): Line 764  class Symbol(object):
764    
765         @param other: object to be added to this object         @param other: object to be added to this object
766         @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}.
767         @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}
768         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
769         """         """
770         return add(self,other)         return add(self,other)
# Line 736  class Symbol(object): Line 775  class Symbol(object):
775    
776         @param other: object this object is added to         @param other: object this object is added to
777         @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}.
778         @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
779         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
780         """         """
781         return add(other,self)         return add(other,self)
# Line 747  class Symbol(object): Line 786  class Symbol(object):
786    
787         @param other: object to be subtracted from this object         @param other: object to be subtracted from this object
788         @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}.
789         @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
790         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
791         """         """
792         return add(self,-other)         return add(self,-other)
# Line 758  class Symbol(object): Line 797  class Symbol(object):
797    
798         @param other: object this object is been subtracted from         @param other: object this object is been subtracted from
799         @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}.
800         @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}.
801         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
802         """         """
803         return add(-self,other)         return add(-self,other)
# Line 769  class Symbol(object): Line 808  class Symbol(object):
808    
809         @param other: object to be mutiplied by this object         @param other: object to be mutiplied by this object
810         @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}.
811         @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}.
812         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
813         """         """
814         return mult(self,other)         return mult(self,other)
# Line 780  class Symbol(object): Line 819  class Symbol(object):
819    
820         @param other: object this object is multiplied with         @param other: object this object is multiplied with
821         @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}.
822         @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.
823         @rtype: L{DependendSymbol} or 0 if other is identical to zero.         @rtype: L{DependendSymbol} or 0 if other is identical to zero.
824         """         """
825         return mult(other,self)         return mult(other,self)
# Line 791  class Symbol(object): Line 830  class Symbol(object):
830    
831         @param other: object dividing this object         @param other: object dividing this object
832         @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}.
833         @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}
834         @rtype: L{DependendSymbol}         @rtype: L{DependendSymbol}
835         """         """
836         return quotient(self,other)         return quotient(self,other)
# Line 802  class Symbol(object): Line 841  class Symbol(object):
841    
842         @param other: object dividing this object         @param other: object dividing this object
843         @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}.
844         @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
845         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
846         """         """
847         return quotient(other,self)         return quotient(other,self)
# Line 813  class Symbol(object): Line 852  class Symbol(object):
852    
853         @param other: exponent         @param other: exponent
854         @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}.
855         @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}
856         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
857         """         """
858         return power(self,other)         return power(self,other)
# Line 824  class Symbol(object): Line 863  class Symbol(object):
863    
864         @param other: basis         @param other: basis
865         @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}.
866         @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
867         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.         @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
868         """         """
869         return power(other,self)         return power(other,self)
870    
871       def __getitem__(self,index):
872           """
873           returns the slice defined by index
874    
875           @param index: defines a
876           @type index: C{slice} or C{int} or a C{tuple} of them
877           @return: a L{Symbol} representing the slice defined by index
878           @rtype: L{DependendSymbol}
879           """
880           return GetSlice_Symbol(self,index)
881    
882  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
883     """     """
884     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.
885     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  
886        
887     Example:     Example::
888        
889     u1=Symbol(shape=(3,4),dim=2,args=[4.])       u1=Symbol(shape=(3,4),dim=2,args=[4.])
890     u2=Symbol(shape=(3,4),dim=2,args=[4.])       u2=Symbol(shape=(3,4),dim=2,args=[4.])
891     print u1==u2       print u1==u2
892     False       False
893        
894        but       but::
895    
896     u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
897     u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
898     u3=DependendSymbol(shape=(2,),dim=2,args=[4.])         u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
899     print u1==u2, u1==u3       print u1==u2, u1==u3
900     True False       True False
901    
902     @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.
903     """     """
# Line 879  class DependendSymbol(Symbol): Line 929  class DependendSymbol(Symbol):
929  #=========================================================  #=========================================================
930  #  Unary operations prserving the shape  #  Unary operations prserving the shape
931  #========================================================  #========================================================
932    class GetSlice_Symbol(DependendSymbol):
933       """
934       L{Symbol} representing getting a slice for a L{Symbol}
935       """
936       def __init__(self,arg,index):
937          """
938          initialization of wherePositive L{Symbol} with argument arg
939          @param arg: argument
940          @type arg: L{Symbol}.
941          @param index: defines index
942          @type index: C{slice} or C{int} or a C{tuple} of them
943          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
944          @raises ValueError: if a step is given
945          """
946          if not isinstance(index,tuple): index=(index,)
947          if len(index)>arg.getRank():
948               raise IndexError,"GetSlice_Symbol: index out of range."
949          sh=()
950          index2=()
951          for i in range(len(index)):
952             ix=index[i]
953             if isinstance(ix,int):
954                if ix<0 or ix>=arg.getShape()[i]:
955                   raise ValueError,"GetSlice_Symbol: index out of range."
956                index2=index2+(ix,)
957             else:
958               if not ix.step==None:
959                 raise ValueError,"GetSlice_Symbol: steping is not supported."
960               if ix.start==None:
961                  s=0
962               else:
963                  s=ix.start
964               if ix.stop==None:
965                  e=arg.getShape()[i]
966               else:
967                  e=ix.stop
968                  if e>arg.getShape()[i]:
969                     raise IndexError,"GetSlice_Symbol: index out of range."
970               index2=index2+(slice(s,e),)
971               if e>s:
972                   sh=sh+(e-s,)
973               elif s>e:
974                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
975          for i in range(len(index),arg.getRank()):
976              index2=index2+(slice(0,arg.getShape()[i]),)
977              sh=sh+(arg.getShape()[i],)
978          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
979    
980       def getMyCode(self,argstrs,format="escript"):
981          """
982          returns a program code that can be used to evaluate the symbol.
983    
984          @param argstrs: gives for each argument a string representing the argument for the evaluation.
985          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
986          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
987          @type format: C{str}
988          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
989          @rtype: C{str}
990          @raise NotImplementedError: if the requested format is not available
991          """
992          if format=="escript" or format=="str"  or format=="text":
993             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
994          else:
995             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
996    
997       def substitute(self,argvals):
998          """
999          assigns new values to symbols in the definition of the symbol.
1000          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1001    
1002          @param argvals: new values assigned to symbols
1003          @type argvals: C{dict} with keywords of type L{Symbol}.
1004          @return: result of the substitution process. Operations are executed as much as possible.
1005          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1006          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1007          """
1008          if argvals.has_key(self):
1009             arg=argvals[self]
1010             if self.isAppropriateValue(arg):
1011                return arg
1012             else:
1013                raise TypeError,"%s: new value is not appropriate."%str(self)
1014          else:
1015             args=self.getSubstitutedArguments(argvals)
1016             arg=args[0]
1017             index=args[1]
1018             return arg.__getitem__(index)
1019    
1020  def log10(arg):  def log10(arg):
1021     """     """
1022     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
1023    
1024     @param arg: argument     @param arg: argument
1025     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1026     @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.
1027     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1028     """     """
1029     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 907  def wherePositive(arg): Line 1045  def wherePositive(arg):
1045    
1046     @param arg: argument     @param arg: argument
1047     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1048     @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.
1049     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1050     """     """
1051     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1052        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1053        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1054        return out        return out
1055     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1056        return arg._wherePositive()        return arg._wherePositive()
# Line 953  class WherePositive_Symbol(DependendSymb Line 1091  class WherePositive_Symbol(DependendSymb
1091        @type format: C{str}        @type format: C{str}
1092        @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.
1093        @rtype: C{str}        @rtype: C{str}
1094        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1095        """        """
1096        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1097            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 989  def whereNegative(arg): Line 1127  def whereNegative(arg):
1127    
1128     @param arg: argument     @param arg: argument
1129     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1130     @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.
1131     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1132     """     """
1133     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1134        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1135        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1136        return out        return out
1137     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1138        return arg._whereNegative()        return arg._whereNegative()
# Line 1035  class WhereNegative_Symbol(DependendSymb Line 1173  class WhereNegative_Symbol(DependendSymb
1173        @type format: C{str}        @type format: C{str}
1174        @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.
1175        @rtype: C{str}        @rtype: C{str}
1176        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1177        """        """
1178        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1179            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1071  def whereNonNegative(arg): Line 1209  def whereNonNegative(arg):
1209    
1210     @param arg: argument     @param arg: argument
1211     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1212     @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.
1213     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1214     """     """
1215     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1216        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1217        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1218        return out        return out
1219     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1220        return arg._whereNonNegative()        return arg._whereNonNegative()
# Line 1101  def whereNonPositive(arg): Line 1239  def whereNonPositive(arg):
1239    
1240     @param arg: argument     @param arg: argument
1241     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1242     @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.
1243     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1244     """     """
1245     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1246        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1247        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1248        return out        return out
1249     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1250        return arg._whereNonPositive()        return arg._whereNonPositive()
# Line 1133  def whereZero(arg,tol=0.): Line 1271  def whereZero(arg,tol=0.):
1271     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1272     @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.
1273     @type tol: C{float}     @type tol: C{float}
1274     @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.
1275     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1276     """     """
1277     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1278        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1279        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1280        return out        return out
1281     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1282        if tol>0.:        return arg._whereZero(tol)
          return whereNegative(abs(arg)-tol)  
       else:  
          return arg._whereZero()  
1283     elif isinstance(arg,float):     elif isinstance(arg,float):
1284        if abs(arg)<=tol:        if abs(arg)<=tol:
1285          return 1.          return 1.
# Line 1182  class WhereZero_Symbol(DependendSymbol): Line 1317  class WhereZero_Symbol(DependendSymbol):
1317        @type format: C{str}        @type format: C{str}
1318        @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.
1319        @rtype: C{str}        @rtype: C{str}
1320        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1321        """        """
1322        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
1323           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])           return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
# Line 1216  def whereNonZero(arg,tol=0.): Line 1351  def whereNonZero(arg,tol=0.):
1351    
1352     @param arg: argument     @param arg: argument
1353     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1354     @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.
1355     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1356     """     """
1357     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1358        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1359        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1360        return out        return out
1361     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1362        if tol>0.:        return arg._whereNonZero(tol)
          return 1.-whereZero(arg,tol)  
       else:  
          return arg._whereNonZero()  
1363     elif isinstance(arg,float):     elif isinstance(arg,float):
1364        if abs(arg)>tol:        if abs(arg)>tol:
1365          return 1.          return 1.
# Line 1243  def whereNonZero(arg,tol=0.): Line 1375  def whereNonZero(arg,tol=0.):
1375     else:     else:
1376        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1377    
1378    def erf(arg):
1379       """
1380       returns erf of argument arg
1381    
1382       @param arg: argument
1383       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1384       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1385       @raises TypeError: if the type of the argument is not expected.
1386       """
1387       if isinstance(arg,escript.Data):
1388          return arg._erf()
1389       else:
1390          raise TypeError,"erf: Unknown argument type."
1391    
1392  def sin(arg):  def sin(arg):
1393     """     """
1394     returns sine of argument arg     returns sine of argument arg
1395    
1396     @param arg: argument     @param arg: argument
1397     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1398     @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.
1399     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1400     """     """
1401     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1287  class Sin_Symbol(DependendSymbol): Line 1433  class Sin_Symbol(DependendSymbol):
1433        @type format: C{str}        @type format: C{str}
1434        @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.
1435        @rtype: C{str}        @rtype: C{str}
1436        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1437        """        """
1438        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1439            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1339  def cos(arg): Line 1485  def cos(arg):
1485    
1486     @param arg: argument     @param arg: argument
1487     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1488     @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.
1489     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1490     """     """
1491     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1377  class Cos_Symbol(DependendSymbol): Line 1523  class Cos_Symbol(DependendSymbol):
1523        @type format: C{str}        @type format: C{str}
1524        @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.
1525        @rtype: C{str}        @rtype: C{str}
1526        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1527        """        """
1528        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1529            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1429  def tan(arg): Line 1575  def tan(arg):
1575    
1576     @param arg: argument     @param arg: argument
1577     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1578     @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.
1579     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1580     """     """
1581     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1467  class Tan_Symbol(DependendSymbol): Line 1613  class Tan_Symbol(DependendSymbol):
1613        @type format: C{str}        @type format: C{str}
1614        @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.
1615        @rtype: C{str}        @rtype: C{str}
1616        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1617        """        """
1618        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1619            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1519  def asin(arg): Line 1665  def asin(arg):
1665    
1666     @param arg: argument     @param arg: argument
1667     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1668     @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.
1669     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1670     """     """
1671     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1557  class Asin_Symbol(DependendSymbol): Line 1703  class Asin_Symbol(DependendSymbol):
1703        @type format: C{str}        @type format: C{str}
1704        @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.
1705        @rtype: C{str}        @rtype: C{str}
1706        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1707        """        """
1708        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1709            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1609  def acos(arg): Line 1755  def acos(arg):
1755    
1756     @param arg: argument     @param arg: argument
1757     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1758     @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.
1759     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1760     """     """
1761     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1647  class Acos_Symbol(DependendSymbol): Line 1793  class Acos_Symbol(DependendSymbol):
1793        @type format: C{str}        @type format: C{str}
1794        @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.
1795        @rtype: C{str}        @rtype: C{str}
1796        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1797        """        """
1798        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1799            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1699  def atan(arg): Line 1845  def atan(arg):
1845    
1846     @param arg: argument     @param arg: argument
1847     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1848     @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.
1849     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1850     """     """
1851     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1737  class Atan_Symbol(DependendSymbol): Line 1883  class Atan_Symbol(DependendSymbol):
1883        @type format: C{str}        @type format: C{str}
1884        @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.
1885        @rtype: C{str}        @rtype: C{str}
1886        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1887        """        """
1888        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1889            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1789  def sinh(arg): Line 1935  def sinh(arg):
1935    
1936     @param arg: argument     @param arg: argument
1937     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1938     @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.
1939     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1940     """     """
1941     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1827  class Sinh_Symbol(DependendSymbol): Line 1973  class Sinh_Symbol(DependendSymbol):
1973        @type format: C{str}        @type format: C{str}
1974        @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.
1975        @rtype: C{str}        @rtype: C{str}
1976        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
1977        """        """
1978        if isinstance(argstrs,list):        if isinstance(argstrs,list):
1979            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1879  def cosh(arg): Line 2025  def cosh(arg):
2025    
2026     @param arg: argument     @param arg: argument
2027     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2028     @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.
2029     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2030     """     """
2031     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 1917  class Cosh_Symbol(DependendSymbol): Line 2063  class Cosh_Symbol(DependendSymbol):
2063        @type format: C{str}        @type format: C{str}
2064        @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.
2065        @rtype: C{str}        @rtype: C{str}
2066        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2067        """        """
2068        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2069            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 1969  def tanh(arg): Line 2115  def tanh(arg):
2115    
2116     @param arg: argument     @param arg: argument
2117     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2118     @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.
2119     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2120     """     """
2121     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2007  class Tanh_Symbol(DependendSymbol): Line 2153  class Tanh_Symbol(DependendSymbol):
2153        @type format: C{str}        @type format: C{str}
2154        @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.
2155        @rtype: C{str}        @rtype: C{str}
2156        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2157        """        """
2158        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2159            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2059  def asinh(arg): Line 2205  def asinh(arg):
2205    
2206     @param arg: argument     @param arg: argument
2207     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2208     @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.
2209     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2210     """     """
2211     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2097  class Asinh_Symbol(DependendSymbol): Line 2243  class Asinh_Symbol(DependendSymbol):
2243        @type format: C{str}        @type format: C{str}
2244        @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.
2245        @rtype: C{str}        @rtype: C{str}
2246        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2247        """        """
2248        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2249            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2149  def acosh(arg): Line 2295  def acosh(arg):
2295    
2296     @param arg: argument     @param arg: argument
2297     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2298     @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.
2299     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2300     """     """
2301     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2187  class Acosh_Symbol(DependendSymbol): Line 2333  class Acosh_Symbol(DependendSymbol):
2333        @type format: C{str}        @type format: C{str}
2334        @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.
2335        @rtype: C{str}        @rtype: C{str}
2336        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2337        """        """
2338        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2339            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2239  def atanh(arg): Line 2385  def atanh(arg):
2385    
2386     @param arg: argument     @param arg: argument
2387     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2388     @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.
2389     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2390     """     """
2391     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2277  class Atanh_Symbol(DependendSymbol): Line 2423  class Atanh_Symbol(DependendSymbol):
2423        @type format: C{str}        @type format: C{str}
2424        @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.
2425        @rtype: C{str}        @rtype: C{str}
2426        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2427        """        """
2428        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2429            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2329  def exp(arg): Line 2475  def exp(arg):
2475    
2476     @param arg: argument     @param arg: argument
2477     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2478     @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.
2479     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2480     """     """
2481     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2367  class Exp_Symbol(DependendSymbol): Line 2513  class Exp_Symbol(DependendSymbol):
2513        @type format: C{str}        @type format: C{str}
2514        @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.
2515        @rtype: C{str}        @rtype: C{str}
2516        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2517        """        """
2518        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2519            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2419  def sqrt(arg): Line 2565  def sqrt(arg):
2565    
2566     @param arg: argument     @param arg: argument
2567     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2568     @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.
2569     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2570     """     """
2571     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2457  class Sqrt_Symbol(DependendSymbol): Line 2603  class Sqrt_Symbol(DependendSymbol):
2603        @type format: C{str}        @type format: C{str}
2604        @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.
2605        @rtype: C{str}        @rtype: C{str}
2606        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2607        """        """
2608        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2609            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2509  def log(arg): Line 2655  def log(arg):
2655    
2656     @param arg: argument     @param arg: argument
2657     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2658     @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.
2659     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2660     """     """
2661     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2547  class Log_Symbol(DependendSymbol): Line 2693  class Log_Symbol(DependendSymbol):
2693        @type format: C{str}        @type format: C{str}
2694        @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.
2695        @rtype: C{str}        @rtype: C{str}
2696        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2697        """        """
2698        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2699            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2599  def sign(arg): Line 2745  def sign(arg):
2745    
2746     @param arg: argument     @param arg: argument
2747     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2748     @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.
2749     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2750     """     """
2751     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2647  class Abs_Symbol(DependendSymbol): Line 2793  class Abs_Symbol(DependendSymbol):
2793        @type format: C{str}        @type format: C{str}
2794        @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.
2795        @rtype: C{str}        @rtype: C{str}
2796        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2797        """        """
2798        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2799            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2699  def minval(arg): Line 2845  def minval(arg):
2845    
2846     @param arg: argument     @param arg: argument
2847     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2848     @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.
2849     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2850     """     """
2851     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2740  class Minval_Symbol(DependendSymbol): Line 2886  class Minval_Symbol(DependendSymbol):
2886        @type format: C{str}        @type format: C{str}
2887        @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.
2888        @rtype: C{str}        @rtype: C{str}
2889        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2890        """        """
2891        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2892            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2776  def maxval(arg): Line 2922  def maxval(arg):
2922    
2923     @param arg: argument     @param arg: argument
2924     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2925     @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.
2926     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
2927     """     """
2928     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
# Line 2817  class Maxval_Symbol(DependendSymbol): Line 2963  class Maxval_Symbol(DependendSymbol):
2963        @type format: C{str}        @type format: C{str}
2964        @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.
2965        @rtype: C{str}        @rtype: C{str}
2966        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
2967        """        """
2968        if isinstance(argstrs,list):        if isinstance(argstrs,list):
2969            argstrs=argstrs[0]            argstrs=argstrs[0]
# Line 2853  def length(arg): Line 2999  def length(arg):
2999    
3000     @param arg: argument     @param arg: argument
3001     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3002     @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.
3003     """     """
3004     return sqrt(inner(arg,arg))     return sqrt(inner(arg,arg))
3005    
# Line 2863  def trace(arg,axis_offset=0): Line 3009  def trace(arg,axis_offset=0):
3009    
3010     @param arg: argument     @param arg: argument
3011     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3012     @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
3013                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3014     @type axis_offset: C{int}     @type axis_offset: C{int}
3015     @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.
3016     @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 2872  def trace(arg,axis_offset=0): Line 3018  def trace(arg,axis_offset=0):
3018     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
3019        sh=arg.shape        sh=arg.shape
3020        if len(sh)<2:        if len(sh)<2:
3021          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3022        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
3023          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
3024        s1=1        s1=1
3025        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
3026        s2=1        s2=1
3027        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3028        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
3029          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)
3030        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))
3031        out=numarray.zeros([s1,s2],numarray.Float)        out=numarray.zeros([s1,s2],numarray.Float64)
3032        for i1 in range(s1):        for i1 in range(s1):
3033          for i2 in range(s2):          for i2 in range(s2):
3034              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
3035        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3036        return out        return out
3037     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3038        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
3039            raise ValueError,"rank of argument must be greater than 1"
3040          if axis_offset<0 or axis_offset>arg.getRank()-2:
3041            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3042          s=list(arg.getShape())        
3043          if not s[axis_offset] == s[axis_offset+1]:
3044            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3045          return arg._trace(axis_offset)
3046     elif isinstance(arg,float):     elif isinstance(arg,float):
3047        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
3048     elif isinstance(arg,int):     elif isinstance(arg,int):
3049        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
3050     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3051        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
3052     else:     else:
3053        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
3054    
 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  
3055  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
3056     """     """
3057     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 2947  class Trace_Symbol(DependendSymbol): Line 3061  class Trace_Symbol(DependendSymbol):
3061        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3062        @param arg: argument of function        @param arg: argument of function
3063        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3064        @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
3065                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3066        @type axis_offset: C{int}        @type axis_offset: C{int}
3067        """        """
3068        if arg.getRank()<2:        if arg.getRank()<2:
3069          raise ValueError,"Trace_Symbol: 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>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3071          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
3072        s=list(arg.getShape())                s=list(arg.getShape())        
3073        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3074          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)
3075        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())
3076    
3077     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 2970  class Trace_Symbol(DependendSymbol): Line 3084  class Trace_Symbol(DependendSymbol):
3084        @type format: C{str}        @type format: C{str}
3085        @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.
3086        @rtype: C{str}        @rtype: C{str}
3087        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3088        """        """
3089        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3090           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])           return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
# Line 3012  class Trace_Symbol(DependendSymbol): Line 3126  class Trace_Symbol(DependendSymbol):
3126        else:        else:
3127           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3128    
3129    def transpose(arg,axis_offset=None):
3130       """
3131       returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3132    
3133       @param arg: argument
3134       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3135       @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.
3136                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3137       @type axis_offset: C{int}
3138       @return: transpose of arg
3139       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3140       """
3141       if isinstance(arg,numarray.NumArray):
3142          if axis_offset==None: axis_offset=int(arg.rank/2)
3143          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3144       elif isinstance(arg,escript.Data):
3145          r=arg.getRank()
3146          if axis_offset==None: axis_offset=int(r/2)
3147          if axis_offset<0 or axis_offset>r:
3148            raise ValueError,"axis_offset must be between 0 and %s"%r
3149          return arg._transpose(axis_offset)
3150       elif isinstance(arg,float):
3151          if not ( axis_offset==0 or axis_offset==None):
3152            raise ValueError,"axis_offset must be 0 for float argument"
3153          return arg
3154       elif isinstance(arg,int):
3155          if not ( axis_offset==0 or axis_offset==None):
3156            raise ValueError,"axis_offset must be 0 for int argument"
3157          return float(arg)
3158       elif isinstance(arg,Symbol):
3159          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3160          return Transpose_Symbol(arg,axis_offset)
3161       else:
3162          raise TypeError,"Unknown argument type."
3163    
3164    class Transpose_Symbol(DependendSymbol):
3165       """
3166       L{Symbol} representing the result of the transpose function
3167       """
3168       def __init__(self,arg,axis_offset=None):
3169          """
3170          initialization of transpose L{Symbol} with argument arg
3171    
3172          @param arg: argument of function
3173          @type arg: L{Symbol}.
3174          @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.
3175                           if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3176          @type axis_offset: C{int}
3177          """
3178          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3179          if axis_offset<0 or axis_offset>arg.getRank():
3180            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3181          s=arg.getShape()
3182          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3183    
3184       def getMyCode(self,argstrs,format="escript"):
3185          """
3186          returns a program code that can be used to evaluate the symbol.
3187    
3188          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3189          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3190          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3191          @type format: C{str}
3192          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3193          @rtype: C{str}
3194          @raise NotImplementedError: if the requested format is not available
3195          """
3196          if format=="escript" or format=="str"  or format=="text":
3197             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3198          else:
3199             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3200    
3201       def substitute(self,argvals):
3202          """
3203          assigns new values to symbols in the definition of the symbol.
3204          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3205    
3206          @param argvals: new values assigned to symbols
3207          @type argvals: C{dict} with keywords of type L{Symbol}.
3208          @return: result of the substitution process. Operations are executed as much as possible.
3209          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3210          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3211          """
3212          if argvals.has_key(self):
3213             arg=argvals[self]
3214             if self.isAppropriateValue(arg):
3215                return arg
3216             else:
3217                raise TypeError,"%s: new value is not appropriate."%str(self)
3218          else:
3219             arg=self.getSubstitutedArguments(argvals)
3220             return transpose(arg[0],axis_offset=arg[1])
3221    
3222       def diff(self,arg):
3223          """
3224          differential of this object
3225    
3226          @param arg: the derivative is calculated with respect to arg
3227          @type arg: L{escript.Symbol}
3228          @return: derivative with respect to C{arg}
3229          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3230          """
3231          if arg==self:
3232             return identity(self.getShape())
3233          else:
3234             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3235    
3236    def swap_axes(arg,axis0=0,axis1=1):
3237       """
3238       returns the swap of arg by swaping the components axis0 and axis1
3239    
3240       @param arg: argument
3241       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3242       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3243       @type axis0: C{int}
3244       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3245       @type axis1: C{int}
3246       @return: C{arg} with swaped components
3247       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3248       """
3249       if axis0 > axis1:
3250          axis0,axis1=axis1,axis0
3251       if isinstance(arg,numarray.NumArray):
3252          return numarray.swapaxes(arg,axis0,axis1)
3253       elif isinstance(arg,escript.Data):
3254          return arg._swap_axes(axis0,axis1)
3255       elif isinstance(arg,float):
3256          raise TyepError,"float argument is not supported."
3257       elif isinstance(arg,int):
3258          raise TyepError,"int argument is not supported."
3259       elif isinstance(arg,Symbol):
3260          return SwapAxes_Symbol(arg,axis0,axis1)
3261       else:
3262          raise TypeError,"Unknown argument type."
3263    
3264    class SwapAxes_Symbol(DependendSymbol):
3265       """
3266       L{Symbol} representing the result of the swap function
3267       """
3268       def __init__(self,arg,axis0=0,axis1=1):
3269          """
3270          initialization of swap L{Symbol} with argument arg
3271    
3272          @param arg: argument
3273          @type arg: L{Symbol}.
3274          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3275          @type axis0: C{int}
3276          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3277          @type axis1: C{int}
3278          """
3279          if arg.getRank()<2:
3280             raise ValueError,"argument must have at least rank 2."
3281          if axis0<0 or axis0>arg.getRank()-1:
3282             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3283          if axis1<0 or axis1>arg.getRank()-1:
3284             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3285          if axis0 == axis1:
3286             raise ValueError,"axis indices must be different."
3287          if axis0 > axis1:
3288             axis0,axis1=axis1,axis0
3289          s=arg.getShape()
3290          s_out=[]
3291          for i in range(len(s)):
3292             if i == axis0:
3293                s_out.append(s[axis1])
3294             elif i == axis1:
3295                s_out.append(s[axis0])
3296             else:
3297                s_out.append(s[i])
3298          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3299    
3300       def getMyCode(self,argstrs,format="escript"):
3301          """
3302          returns a program code that can be used to evaluate the symbol.
3303    
3304          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3305          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3306          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3307          @type format: C{str}
3308          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3309          @rtype: C{str}
3310          @raise NotImplementedError: if the requested format is not available
3311          """
3312          if format=="escript" or format=="str"  or format=="text":
3313             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3314          else:
3315             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3316    
3317       def substitute(self,argvals):
3318          """
3319          assigns new values to symbols in the definition of the symbol.
3320          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3321    
3322          @param argvals: new values assigned to symbols
3323          @type argvals: C{dict} with keywords of type L{Symbol}.
3324          @return: result of the substitution process. Operations are executed as much as possible.
3325          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3326          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3327          """
3328          if argvals.has_key(self):
3329             arg=argvals[self]
3330             if self.isAppropriateValue(arg):
3331                return arg
3332             else:
3333                raise TypeError,"%s: new value is not appropriate."%str(self)
3334          else:
3335             arg=self.getSubstitutedArguments(argvals)
3336             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3337    
3338       def diff(self,arg):
3339          """
3340          differential of this object
3341    
3342          @param arg: the derivative is calculated with respect to arg
3343          @type arg: L{escript.Symbol}
3344          @return: derivative with respect to C{arg}
3345          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3346          """
3347          if arg==self:
3348             return identity(self.getShape())
3349          else:
3350             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3351    
3352    def symmetric(arg):
3353        """
3354        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3355    
3356        @param arg: square matrix. Must have rank 2 or 4 and be square.
3357        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3358        @return: symmetric part of arg
3359        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3360        """
3361        if isinstance(arg,numarray.NumArray):
3362          if arg.rank==2:
3363            if not (arg.shape[0]==arg.shape[1]):
3364               raise ValueError,"argument must be square."
3365          elif arg.rank==4:
3366            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3367               raise ValueError,"argument must be square."
3368          else:
3369            raise ValueError,"rank 2 or 4 is required."
3370          return (arg+transpose(arg))/2
3371        elif isinstance(arg,escript.Data):
3372          if arg.getRank()==2:
3373            if not (arg.getShape()[0]==arg.getShape()[1]):
3374               raise ValueError,"argument must be square."
3375            return arg._symmetric()
3376          elif arg.getRank()==4:
3377            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3378               raise ValueError,"argument must be square."
3379            return arg._symmetric()
3380          else:
3381            raise ValueError,"rank 2 or 4 is required."
3382        elif isinstance(arg,float):
3383          return arg
3384        elif isinstance(arg,int):
3385          return float(arg)
3386        elif isinstance(arg,Symbol):
3387          if arg.getRank()==2:
3388            if not (arg.getShape()[0]==arg.getShape()[1]):
3389               raise ValueError,"argument must be square."
3390          elif arg.getRank()==4:
3391            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3392               raise ValueError,"argument must be square."
3393          else:
3394            raise ValueError,"rank 2 or 4 is required."
3395          return (arg+transpose(arg))/2
3396        else:
3397          raise TypeError,"symmetric: Unknown argument type."
3398    
3399    def nonsymmetric(arg):
3400        """
3401        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3402    
3403        @param arg: square matrix. Must have rank 2 or 4 and be square.
3404        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3405        @return: nonsymmetric part of arg
3406        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3407        """
3408        if isinstance(arg,numarray.NumArray):
3409          if arg.rank==2:
3410            if not (arg.shape[0]==arg.shape[1]):
3411               raise ValueError,"nonsymmetric: argument must be square."
3412          elif arg.rank==4:
3413            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3414               raise ValueError,"nonsymmetric: argument must be square."
3415          else:
3416            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3417          return (arg-transpose(arg))/2
3418        elif isinstance(arg,escript.Data):
3419          if arg.getRank()==2:
3420            if not (arg.getShape()[0]==arg.getShape()[1]):
3421               raise ValueError,"argument must be square."
3422            return arg._nonsymmetric()
3423          elif arg.getRank()==4:
3424            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3425               raise ValueError,"argument must be square."
3426            return arg._nonsymmetric()
3427          else:
3428            raise ValueError,"rank 2 or 4 is required."
3429        elif isinstance(arg,float):
3430          return arg
3431        elif isinstance(arg,int):
3432          return float(arg)
3433        elif isinstance(arg,Symbol):
3434          if arg.getRank()==2:
3435            if not (arg.getShape()[0]==arg.getShape()[1]):
3436               raise ValueError,"nonsymmetric: argument must be square."
3437          elif arg.getRank()==4:
3438            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3439               raise ValueError,"nonsymmetric: argument must be square."
3440          else:
3441            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3442          return (arg-transpose(arg))/2
3443        else:
3444          raise TypeError,"nonsymmetric: Unknown argument type."
3445    
3446  def inverse(arg):  def inverse(arg):
3447      """      """
3448      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3449    
3450      @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.
3451      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3452      @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])
3453      @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
3454        @note: for L{escript.Data} objects the dimension is restricted to 3.
3455      """      """
3456        import numarray.linear_algebra # This statement should be after the next statement but then somehow numarray is gone.
3457      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3458        return numarray.linear_algebra.inverse(arg)        return numarray.linear_algebra.inverse(arg)
3459      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
# Line 3113  class Inverse_Symbol(DependendSymbol): Line 3546  class Inverse_Symbol(DependendSymbol):
3546        @type format: C{str}        @type format: C{str}
3547        @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.
3548        @rtype: C{str}        @rtype: C{str}
3549        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3550        """        """
3551        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
3552           return "inverse(%s)"%argstrs[0]           return "inverse(%s)"%argstrs[0]
# Line 3153  class Inverse_Symbol(DependendSymbol): Line 3586  class Inverse_Symbol(DependendSymbol):
3586        if arg==self:        if arg==self:
3587           return identity(self.getShape())           return identity(self.getShape())
3588        else:        else:
3589           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3590    
3591    def eigenvalues(arg):
3592        """
3593        returns the eigenvalues of the square matrix arg.
3594    
3595        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3596                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3597        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3598        @return: the eigenvalues in increasing order.
3599        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3600        @note: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3601        """
3602        if isinstance(arg,numarray.NumArray):
3603          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3604          out.sort()
3605          return out
3606        elif isinstance(arg,escript.Data):
3607          return arg._eigenvalues()
3608        elif isinstance(arg,Symbol):
3609          if not arg.getRank()==2:
3610            raise ValueError,"eigenvalues: argument must have rank 2"
3611          s=arg.getShape()      
3612          if not s[0] == s[1]:
3613            raise ValueError,"eigenvalues: argument must be a square matrix."
3614          if s[0]==1:
3615              return arg[0]
3616          elif s[0]==2:
3617              arg1=symmetric(arg)
3618              A11=arg1[0,0]
3619              A12=arg1[0,1]
3620              A22=arg1[1,1]
3621              trA=(A11+A22)/2.
3622              A11-=trA
3623              A22-=trA
3624              s=sqrt(A12**2-A11*A22)
3625              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3626          elif s[0]==3:
3627              arg1=symmetric(arg)
3628              A11=arg1[0,0]
3629              A12=arg1[0,1]
3630              A22=arg1[1,1]
3631              A13=arg1[0,2]
3632              A23=arg1[1,2]
3633              A33=arg1[2,2]
3634              trA=(A11+A22+A33)/3.
3635              A11-=trA
3636              A22-=trA
3637              A33-=trA
3638              A13_2=A13**2
3639              A23_2=A23**2
3640              A12_2=A12**2
3641              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3642              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3643              sq_p=sqrt(p/3.)
3644              alpha_3=acos(clip(-q*(sq_p+whereZero(p,0.)*1.e-15)**(-3.)/2.,-1.,1.))/3.  # whereZero is protection against divison by zero
3645              sq_p*=2.
3646              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3647               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3648               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3649              return trA+sq_p*f
3650          else:
3651             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3652        elif isinstance(arg,float):
3653          return arg
3654        elif isinstance(arg,int):
3655          return float(arg)
3656        else:
3657          raise TypeError,"eigenvalues: Unknown argument type."
3658    
3659    def eigenvalues_and_eigenvectors(arg):
3660        """
3661        returns the eigenvalues and eigenvectors of the square matrix arg.
3662    
3663        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3664                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3665        @type arg: L{escript.Data}
3666        @return: the eigenvalues and eigenvectors. The eigenvalues are ordered by increasing value. The
3667                 eigenvectors are orthogonal and normalized. If V are the eigenvectors than V[:,i] is
3668                 the eigenvector coresponding to the i-th eigenvalue.
3669        @rtype: L{tuple} of L{escript.Data}.
3670        @note: The dimension is restricted to 3.
3671        """
3672        if isinstance(arg,numarray.NumArray):
3673          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3674        elif isinstance(arg,escript.Data):
3675          return arg._eigenvalues_and_eigenvectors()
3676        elif isinstance(arg,Symbol):
3677          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3678        elif isinstance(arg,float):
3679          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3680        elif isinstance(arg,int):
3681          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3682        else:
3683          raise TypeError,"eigenvalues: Unknown argument type."
3684  #=======================================================  #=======================================================
3685  #  Binary operations:  #  Binary operations:
3686  #=======================================================  #=======================================================
# Line 3197  class Add_Symbol(DependendSymbol): Line 3724  class Add_Symbol(DependendSymbol):
3724         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3725         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3726         """         """
3727         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3728         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3729         if not sh0==sh1:         if not sh0==sh1:
3730            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3731         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 3213  class Add_Symbol(DependendSymbol): Line 3740  class Add_Symbol(DependendSymbol):
3740        @type format: C{str}        @type format: C{str}
3741        @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.
3742        @rtype: C{str}        @rtype: C{str}
3743        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3744        """        """
3745        if format=="str" or format=="text":        if format=="str" or format=="text":
3746           return "(%s)+(%s)"%(argstrs[0],argstrs[1])           return "(%s)+(%s)"%(argstrs[0],argstrs[1])
# Line 3272  def mult(arg0,arg1): Line 3799  def mult(arg0,arg1):
3799         """         """
3800         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3801         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3802            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3803         else:         else:
3804            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3805                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3296  class Mult_Symbol(DependendSymbol): Line 3823  class Mult_Symbol(DependendSymbol):
3823         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3824         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3825         """         """
3826         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3827         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3828         if not sh0==sh1:         if not sh0==sh1:
3829            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3830         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 3312  class Mult_Symbol(DependendSymbol): Line 3839  class Mult_Symbol(DependendSymbol):
3839        @type format: C{str}        @type format: C{str}
3840        @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.
3841        @rtype: C{str}        @rtype: C{str}
3842        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3843        """        """
3844        if format=="str" or format=="text":        if format=="str" or format=="text":
3845           return "(%s)*(%s)"%(argstrs[0],argstrs[1])           return "(%s)*(%s)"%(argstrs[0],argstrs[1])
# Line 3372  def quotient(arg0,arg1): Line 3899  def quotient(arg0,arg1):
3899         """         """
3900         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3901         if testForZero(args[0]):         if testForZero(args[0]):
3902            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3903         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3904            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3905               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3401  class Quotient_Symbol(DependendSymbol): Line 3928  class Quotient_Symbol(DependendSymbol):
3928         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3929         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3930         """         """
3931         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3932         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3933         if not sh0==sh1:         if not sh0==sh1:
3934            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3935         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 3417  class Quotient_Symbol(DependendSymbol): Line 3944  class Quotient_Symbol(DependendSymbol):
3944        @type format: C{str}        @type format: C{str}
3945        @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.
3946        @rtype: C{str}        @rtype: C{str}
3947        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
3948        """        """
3949        if format=="str" or format=="text":        if format=="str" or format=="text":
3950           return "(%s)/(%s)"%(argstrs[0],argstrs[1])           return "(%s)/(%s)"%(argstrs[0],argstrs[1])
# Line 3478  def power(arg0,arg1): Line 4005  def power(arg0,arg1):
4005         """         """
4006         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4007         if testForZero(args[0]):         if testForZero(args[0]):
4008            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4009         elif testForZero(args[1]):         elif testForZero(args[1]):
4010            return numarray.ones(args[0],numarray.Float)            return numarray.ones(getShape(args[1]),numarray.Float64)
4011         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4012            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4013         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 3503  class Power_Symbol(DependendSymbol): Line 4030  class Power_Symbol(DependendSymbol):
4030         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4031         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4032         """         """
4033         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4034         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4035         if not sh0==sh1:         if not sh0==sh1:
4036            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4037         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3521  class Power_Symbol(DependendSymbol): Line 4048  class Power_Symbol(DependendSymbol):
4048        @type format: C{str}        @type format: C{str}
4049        @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.
4050        @rtype: C{str}        @rtype: C{str}
4051        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4052        """        """
4053        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4054           return "(%s)**(%s)"%(argstrs[0],argstrs[1])           return "(%s)**(%s)"%(argstrs[0],argstrs[1])
# Line 3603  def minimum(*args): Line 4130  def minimum(*args):
4130            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4131      return out      return out
4132    
4133  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4134      """      """
4135      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4136    
4137      @param arg: argument      @param arg: argument
4138      @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}
4139      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4140      @type arg: C{float}      @type minval: C{float} or C{None}
4141      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4142      @type arg: C{float}      @type maxval: C{float} or C{None}
4143      @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
4144               less then maxval are unchanged.               less then maxval are unchanged.
4145      @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
4146      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4147      """      """
4148      if minval>maxval:      if not minval==None and not maxval==None:
4149         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4150      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4151        if minval == None:
4152            tmp=arg
4153        else:
4154            tmp=maximum(minval,arg)
4155        if maxval == None:
4156            return tmp
4157        else:
4158            return minimum(tmp,maxval)
4159    
4160        
4161  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3637  def inner(arg0,arg1): Line 4172  def inner(arg0,arg1):
4172      @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}
4173      @param arg1: second argument      @param arg1: second argument
4174      @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}
4175      @return : the inner product of arg0 and arg1 at each data point      @return: the inner product of arg0 and arg1 at each data point
4176      @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
4177      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4178      """      """
4179      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4180      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4181      if not sh0==sh1:      if not sh0==sh1:
4182          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4183      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4184    
4185    def outer(arg0,arg1):
4186        """
4187        the outer product of the two argument:
4188    
4189        out[t,s]=arg0[t]*arg1[s]
4190    
4191        where
4192    
4193            - s runs through arg0.Shape
4194            - t runs through arg1.Shape
4195    
4196        @param arg0: first argument
4197        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4198        @param arg1: second argument
4199        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4200        @return: the outer product of arg0 and arg1 at each data point
4201        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4202        """
4203        return generalTensorProduct(arg0,arg1,axis_offset=0)
4204    
4205  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4206      """      """
4207        see L{matrix_mult}
4208        """
4209        return matrix_mult(arg0,arg1)
4210    
4211    def matrix_mult(arg0,arg1):
4212        """
4213      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4214    
4215      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
4216    
4217            or      or
4218    
4219      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4220    
4221      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.
4222    
4223      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4224      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3667  def matrixmult(arg0,arg1): Line 4228  def matrixmult(arg0,arg1):
4228      @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
4229      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4230      """      """
4231      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4232      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4233      if not len(sh0)==2 :      if not len(sh0)==2 :
4234          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4235      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4236          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4237      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4238    
4239  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4240      """      """
4241      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  
4242      """      """
4243      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
   
4244    
4245  def tensormult(arg0,arg1):  def tensor_mult(arg0,arg1):
4246      """      """
4247      the tensor product of the two argument:      the tensor product of the two argument:
   
4248            
4249      for arg0 of rank 2 this is      for arg0 of rank 2 this is
4250    
4251      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]        out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]  
4252    
4253                   or      or
4254    
4255      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4256    
# Line 3712  def tensormult(arg0,arg1): Line 4259  def tensormult(arg0,arg1):
4259    
4260      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]
4261                                
4262                   or      or
4263    
4264      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]
4265    
4266                   or      or
4267    
4268      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]
4269    
4270      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  
4271      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.
4272    
4273      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4274      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 3730  def tensormult(arg0,arg1): Line 4277  def tensormult(arg0,arg1):
4277      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4278      @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
4279      """      """
4280      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4281      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4282      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4283         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4284      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):
4285         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4286      else:      else:
4287          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4288    
4289  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4290      """      """
# Line 3745  def generalTensorProduct(arg0,arg1,axis_ Line 4292  def generalTensorProduct(arg0,arg1,axis_
4292    
4293      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]      out[s,t]=S{Sigma}_r arg0[s,r]*arg1[r,t]
4294    
4295      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:]  
4296    
4297      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]
4298      in the second case the two last dimensions of arg0 must match the shape of arg1.          - r runs trough arg0.Shape[:axis_offset]
4299            - t runs through arg1.Shape[axis_offset:]
4300    
4301      @param arg0: first argument      @param arg0: first argument
4302      @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}
4303      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4304      @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}
4305      @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.
4306      @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 3767  def generalTensorProduct(arg0,arg1,axis_ Line 4313  def generalTensorProduct(arg0,arg1,axis_
4313             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4314         else:         else:
4315             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4316                 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)
4317             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4318             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4319             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 3777  def generalTensorProduct(arg0,arg1,axis_ Line 4323  def generalTensorProduct(arg0,arg1,axis_
4323             for i in sh1[:axis_offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4324             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4325             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4326             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4327             for i0 in range(d0):             for i0 in range(d0):
4328                      for i1 in range(d1):                      for i1 in range(d1):
4329                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
# Line 3793  def generalTensorProduct(arg0,arg1,axis_ Line 4339  def generalTensorProduct(arg0,arg1,axis_
4339                                    
4340  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4341     """     """
4342     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4343     """     """
4344     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4345         """         """
4346         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4347    
4348         @param arg0: numerator         @param arg0: first argument
4349         @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}.
4350         @param arg1: denominator         @param arg1: second argument
4351         @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}.
4352         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4353         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4354         """         """
4355         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4356         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4357         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4358         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4359         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 3826  class GeneralTensorProduct_Symbol(Depend Line 4372  class GeneralTensorProduct_Symbol(Depend
4372        @type format: C{str}        @type format: C{str}
4373        @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.
4374        @rtype: C{str}        @rtype: C{str}
4375        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4376        """        """
4377        if format=="escript" or format=="str" or format=="text":        if format=="escript" or format=="str" or format=="text":
4378           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 3854  class GeneralTensorProduct_Symbol(Depend Line 4400  class GeneralTensorProduct_Symbol(Depend
4400           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4401           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4402    
4403  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4404      "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!!!"
4405      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4406      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
4407      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  def transposed_matrix_mult(arg0,arg1):
4408      shape10=arg1.getShape()[:axis_offset]      """
4409      shape1=arg1.getShape()[axis_offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4410      if not shape01==shape10:  
4411          raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
   
     # whatr function space should be used? (this here is not good!)  
     fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
     # create return value:  
     out=escript.Data(0.,tuple(shape0+shape1),fs)  
     #  
     s0=[[]]  
     for k in shape0:  
           s=[]  
           for j in s0:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s0=s  
     s1=[[]]  
     for k in shape1:  
           s=[]  
           for j in s1:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s1=s  
     s01=[[]]  
     for k in shape01:  
           s=[]  
           for j in s01:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s01=s  
   
     for i0 in s0:  
        for i1 in s1:  
          s=escript.Scalar(0.,fs)  
          for i01 in s01:  
             s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))  
          out.__setitem__(tuple(i0+i1),s)  
     return out  
4412    
4413        or
4414    
4415        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4416    
4417        The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4418    
4419        The first dimension of arg0 and arg1 must match.
4420    
4421        @param arg0: first argument of rank 2
4422        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4423        @param arg1: second argument of at least rank 1
4424        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4425        @return: the product of the transposed of arg0 and arg1 at each data point
4426        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4427        @raise ValueError: if the shapes of the arguments are not appropriate
4428        """
4429        sh0=getShape(arg0)
4430        sh1=getShape(arg1)
4431        if not len(sh0)==2 :
4432            raise ValueError,"first argument must have rank 2"
4433        if not len(sh1)==2 and not len(sh1)==1:
4434            raise ValueError,"second argument must have rank 1 or 2"
4435        return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4436    
4437    def transposed_tensor_mult(arg0,arg1):
4438        """
4439        the tensor product of the transposed of the first and the second argument
4440        
4441        for arg0 of rank 2 this is
4442    
4443        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4444    
4445        or
4446    
4447        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4448    
4449      
4450        and for arg0 of rank 4 this is
4451    
4452        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4453                  
4454        or
4455    
4456        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4457    
4458        or
4459    
4460        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4461    
4462        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4463        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4464    
4465        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4466    
4467        @param arg0: first argument of rank 2 or 4
4468        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4469        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4470        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4471        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4472        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4473        """
4474        sh0=getShape(arg0)
4475        sh1=getShape(arg1)
4476        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4477           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4478        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4479           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4480        else:
4481            raise ValueError,"first argument must have rank 2 or 4"
4482    
4483    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4484        """
4485        generalized tensor product of transposed of arg0 and arg1:
4486    
4487        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4488    
4489        where
4490    
4491            - s runs through arg0.Shape[axis_offset:]
4492            - r runs trough arg0.Shape[:axis_offset]
4493            - t runs through arg1.Shape[axis_offset:]
4494    
4495        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4496        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4497    
4498        @param arg0: first argument
4499        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4500        @param arg1: second argument
4501        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4502        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4503        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4504        """
4505        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4506        arg0,arg1=matchType(arg0,arg1)
4507        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4508        if isinstance(arg0,numarray.NumArray):
4509           if isinstance(arg1,Symbol):
4510               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4511           else:
4512               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4513                   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)
4514               arg0_c=arg0.copy()
4515               arg1_c=arg1.copy()
4516               sh0,sh1=arg0.shape,arg1.shape
4517               d0,d1,d01=1,1,1
4518               for i in sh0[axis_offset:]: d0*=i
4519               for i in sh1[axis_offset:]: d1*=i
4520               for i in sh0[:axis_offset]: d01*=i
4521               arg0_c.resize((d01,d0))
4522               arg1_c.resize((d01,d1))
4523               out=numarray.zeros((d0,d1),numarray.Float64)
4524               for i0 in range(d0):
4525                        for i1 in range(d1):
4526                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4527               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4528               return out
4529        elif isinstance(arg0,escript.Data):
4530           if isinstance(arg1,Symbol):
4531               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4532           else:
4533               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4534        else:      
4535           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4536                    
4537    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4538       """
4539       Symbol representing the general tensor product of the transposed of arg0 and arg1
4540       """
4541       def __init__(self,arg0,arg1,axis_offset=0):
4542           """
4543           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4544    
4545           @param arg0: first argument
4546           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4547           @param arg1: second argument
4548           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4549           @raise ValueError: inconsistent dimensions of arguments.
4550           @note: if both arguments have a spatial dimension, they must equal.
4551           """
4552           sh_arg0=getShape(arg0)
4553           sh_arg1=getShape(arg1)
4554           sh01=sh_arg0[:axis_offset]
4555           sh10=sh_arg1[:axis_offset]
4556           sh0=sh_arg0[axis_offset:]
4557           sh1=sh_arg1[axis_offset:]
4558           if not sh01==sh10:
4559               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)
4560           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4561    
4562       def getMyCode(self,argstrs,format="escript"):
4563          """
4564          returns a program code that can be used to evaluate the symbol.
4565    
4566          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4567          @type argstrs: C{list} of length 2 of C{str}.
4568          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4569          @type format: C{str}
4570          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4571          @rtype: C{str}
4572          @raise NotImplementedError: if the requested format is not available
4573          """
4574          if format=="escript" or format=="str" or format=="text":
4575             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4576          else:
4577             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4578    
4579       def substitute(self,argvals):
4580          """
4581          assigns new values to symbols in the definition of the symbol.
4582          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4583    
4584          @param argvals: new values assigned to symbols
4585          @type argvals: C{dict} with keywords of type L{Symbol}.
4586          @return: result of the substitution process. Operations are executed as much as possible.
4587          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4588          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4589          """
4590          if argvals.has_key(self):
4591             arg=argvals[self]
4592             if self.isAppropriateValue(arg):
4593                return arg
4594             else:
4595                raise TypeError,"%s: new value is not appropriate."%str(self)
4596          else:
4597             args=self.getSubstitutedArguments(argvals)
4598             return generalTransposedTensorProduct(args[0],args[1],args[2])
4599    
4600    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4601        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4602        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4603    
4604    def matrix_transposed_mult(arg0,arg1):
4605        """
4606        matrix-transposed(matrix) product of the two argument:
4607    
4608        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4609    
4610        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4611    
4612        The last dimensions of arg0 and arg1 must match.
4613    
4614        @param arg0: first argument of rank 2
4615        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4616        @param arg1: second argument of rank 2
4617        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4618        @return: the product of arg0 and the transposed of arg1 at each data point
4619        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4620        @raise ValueError: if the shapes of the arguments are not appropriate
4621        """
4622        sh0=getShape(arg0)
4623        sh1=getShape(arg1)
4624        if not len(sh0)==2 :
4625            raise ValueError,"first argument must have rank 2"
4626        if not len(sh1)==2 and not len(sh1)==1:
4627            raise ValueError,"second argument must have rank 1 or 2"
4628        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4629    
4630    def tensor_transposed_mult(arg0,arg1):
4631        """
4632        the tensor product of the first and the transpose of the second argument
4633        
4634        for arg0 of rank 2 this is
4635    
4636        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4637    
4638        and for arg0 of rank 4 this is
4639    
4640        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4641                  
4642        or
4643    
4644        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4645    
4646        In the first case the the second dimension of arg0 and arg1 must match and  
4647        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4648    
4649        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4650    
4651        @param arg0: first argument of rank 2 or 4
4652        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4653        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4654        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4655        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4656        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4657        """
4658        sh0=getShape(arg0)
4659        sh1=getShape(arg1)
4660        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4661           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4662        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4663           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4664        else:
4665            raise ValueError,"first argument must have rank 2 or 4"
4666    
4667    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4668        """
4669        generalized tensor product of transposed of arg0 and arg1:
4670    
4671        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4672    
4673        where
4674    
4675            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4676            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4677            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4678    
4679        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4680        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4681    
4682        @param arg0: first argument
4683        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4684        @param arg1: second argument
4685        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4686        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4687        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4688        """
4689        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4690        arg0,arg1=matchType(arg0,arg1)
4691        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4692        if isinstance(arg0,numarray.NumArray):
4693           if isinstance(arg1,Symbol):
4694               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4695           else:
4696               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4697                   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)
4698               arg0_c=arg0.copy()
4699               arg1_c=arg1.copy()
4700               sh0,sh1=arg0.shape,arg1.shape
4701               d0,d1,d01=1,1,1
4702               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4703               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4704               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4705               arg0_c.resize((d0,d01))
4706               arg1_c.resize((d1,d01))
4707               out=numarray.zeros((d0,d1),numarray.Float64)
4708               for i0 in range(d0):
4709                        for i1 in range(d1):
4710                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4711               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4712               return out
4713        elif isinstance(arg0,escript.Data):
4714           if isinstance(arg1,Symbol):
4715               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4716           else:
4717               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4718        else:      
4719           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4720                    
4721    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4722       """
4723       Symbol representing the general tensor product of arg0 and the transpose of arg1
4724       """
4725       def __init__(self,arg0,arg1,axis_offset=0):
4726           """
4727           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4728    
4729           @param arg0: first argument
4730           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4731           @param arg1: second argument
4732           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4733           @raise ValueError: inconsistent dimensions of arguments.
4734           @note: if both arguments have a spatial dimension, they must equal.
4735           """
4736           sh_arg0=getShape(arg0)
4737           sh_arg1=getShape(arg1)
4738           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4739           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4740           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4741           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4742           if not sh01==sh10:
4743               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)
4744           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4745    
4746       def getMyCode(self,argstrs,format="escript"):
4747          """
4748          returns a program code that can be used to evaluate the symbol.
4749    
4750          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4751          @type argstrs: C{list} of length 2 of C{str}.
4752          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4753          @type format: C{str}
4754          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4755          @rtype: C{str}
4756          @raise NotImplementedError: if the requested format is not available
4757          """
4758          if format=="escript" or format=="str" or format=="text":
4759             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4760          else:
4761             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4762    
4763       def substitute(self,argvals):
4764          """
4765          assigns new values to symbols in the definition of the symbol.
4766          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4767    
4768          @param argvals: new values assigned to symbols
4769          @type argvals: C{dict} with keywords of type L{Symbol}.
4770          @return: result of the substitution process. Operations are executed as much as possible.
4771          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4772          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4773          """
4774          if argvals.has_key(self):
4775             arg=argvals[self]
4776             if self.isAppropriateValue(arg):
4777                return arg
4778             else:
4779                raise TypeError,"%s: new value is not appropriate."%str(self)
4780          else:
4781             args=self.getSubstitutedArguments(argvals)
4782             return generalTensorTransposedProduct(args[0],args[1],args[2])
4783    
4784    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4785        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4786        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4787    
4788  #=========================================================  #=========================================================
4789  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
# Line 3917  def grad(arg,where=None): Line 4805  def grad(arg,where=None):
4805                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4806      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4807      @return: gradient of arg.      @return: gradient of arg.
4808      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
4809      """      """
4810      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4811         return Grad_Symbol(arg,where)         return Grad_Symbol(arg,where)
# Line 3957  class Grad_Symbol(DependendSymbol): Line 4845  class Grad_Symbol(DependendSymbol):
4845        @type format: C{str}        @type format: C{str}
4846        @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.
4847        @rtype: C{str}        @rtype: C{str}
4848        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4849        """        """
4850        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4851           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])           return "grad(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4010  def integrate(arg,where=None): Line 4898  def integrate(arg,where=None):
4898                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
4899      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
4900      @return: integral of arg.      @return: integral of arg.
4901      @rtype:  C{float}, C{numarray.NumArray} or L{Symbol}      @rtype: C{float}, C{numarray.NumArray} or L{Symbol}
4902      """      """
4903      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4904         return Integrate_Symbol(arg,where)         return Integrate_Symbol(arg,where)
# Line 4048  class Integrate_Symbol(DependendSymbol): Line 4936  class Integrate_Symbol(DependendSymbol):
4936        @type format: C{str}        @type format: C{str}
4937        @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.
4938        @rtype: C{str}        @rtype: C{str}
4939        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
4940        """        """
4941        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
4942           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "integrate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4100  def interpolate(arg,where): Line 4988  def interpolate(arg,where):
4988      @param where: FunctionSpace to be interpolated to      @param where: FunctionSpace to be interpolated to
4989      @type where: L{escript.FunctionSpace}      @type where: L{escript.FunctionSpace}
4990      @return: interpolated argument      @return: interpolated argument
4991      @rtype:  C{escript.Data} or L{Symbol}      @rtype: C{escript.Data} or L{Symbol}
4992      """      """
4993      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
4994         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
# Line 4131  class Interpolate_Symbol(DependendSymbol Line 5019  class Interpolate_Symbol(DependendSymbol
5019        @type format: C{str}        @type format: C{str}
5020        @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.
5021        @rtype: C{str}        @rtype: C{str}
5022        @raise: NotImplementedError: if the requested format is not available        @raise NotImplementedError: if the requested format is not available
5023        """        """
5024        if format=="escript" or format=="str"  or format=="text":        if format=="escript" or format=="str"  or format=="text":
5025           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])           return "interpolate(%s,where=%s)"%(argstrs[0],argstrs[1])
# Line 4184  def div(arg,where=None): Line 5072  def div(arg,where=None):
5072                    If not present or C{None} an appropriate default is used.                    If not present or C{None} an appropriate default is used.
5073      @type where: C{None} or L{escript.FunctionSpace}      @type where: C{None} or L{escript.FunctionSpace}
5074      @return: divergence of arg.      @return: divergence of arg.
5075      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5076      """      """
5077      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5078          dim=arg.getDim()          dim=arg.getDim()
# Line 4206  def jump(arg,domain=None): Line 5094  def jump(arg,domain=None):
5094                     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.
5095      @type domain: C{None} or L{escript.Domain}      @type domain: C{None} or L{escript.Domain}
5096      @return: jump of arg      @return: jump of arg
5097      @rtype:  L{escript.Data} or L{Symbol}      @rtype: L{escript.Data} or L{Symbol}
5098      """      """
5099      if domain==None: domain=arg.getDomain()      if domain==None: domain=arg.getDomain()
5100      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))      return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
# Line 4218  def L2(arg): Line 5106  def L2(arg):
5106      @param arg: function which L2 to be calculated.      @param arg: function which L2 to be calculated.
5107      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
5108      @return: L2 norm of arg.      @return: L2 norm of arg.
5109      @rtype:  L{float} or L{Symbol}      @rtype: L{float} or L{Symbol}
5110      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5111      """      """
5112      return sqrt(integrate(inner(arg,arg)))      return sqrt(integrate(inner(arg,arg)))
5113  #=============================  #=============================
5114  #  #
 # wrapper for various functions: if the argument has attribute the function name  
 # as an argument it calls the corresponding methods. Otherwise the corresponding  
 # numarray function is called.  
   
 # functions involving the underlying Domain:  
   
 def transpose(arg,axis=None):  
     """  
     Returns the transpose of the Data object arg.  
   
     @param arg:  
     """  
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
     if isinstance(arg,Symbol):  
        return Transpose_Symbol(arg,axis=r)  
     if isinstance(arg,escript.Data):  
        # hack for transpose  
        r=arg.getRank()  
        if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"  
        s=arg.getShape()  
        out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())  
        for i in range(s[0]):  
           for j in range(s[1]):  
              out[j,i]=arg[i,j]  
        return out  
        # end hack for transpose  
        return arg.transpose(axis)  
     else:  
        return numarray.transpose(arg,axis=axis)  
   
   
5115    
5116  def reorderComponents(arg,index):  def reorderComponents(arg,index):
5117      """      """
5118      resorts the component of arg according to index      resorts the component of arg according to index
5119    
5120      """      """
5121      pass      raise NotImplementedError
5122  #  #
5123  # $Log: util.py,v $  # $Log: util.py,v $
5124  # Revision 1.14.2.16  2005/10/19 06:09:57  gross  # Revision 1.14.2.16  2005/10/19 06:09:57  gross

Legend:
Removed from v.443  
changed lines
  Added in v.1247

  ViewVC Help
Powered by ViewVC 1.1.26