/[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 795 by ksteube, Mon Jul 31 01:23:58 2006 UTC revision 1809 by ksteube, Thu Sep 25 06:43:44 2008 UTC
# Line 1  Line 1 
1  # $Id$  
2    ########################################################
3    #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
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.uq.edu.au/esscc/escript-finley"
21    
22  """  """
23  Utility functions for escript  Utility functions for escript
# Line 9  Utility functions for escript Line 28  Utility functions for escript
28  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
29  @var __version__: version  @var __version__: version
30  @var __date__: date of the version  @var __date__: date of the version
31    @var EPSILON: smallest positive value with 1.<1.+EPSILON
32  """  """
33                                                                                                                                                                                                        
34  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys/escript"  
 __version__="$Revision$"  
 __date__="$Date$"  
35    
36    
37  import math  import math
38  import numarray  import numarray
39  import escript  import escript
40  import os  import os
41    from esys.escript import C_GeneralTensorProduct
42    from esys.escript import getVersion
43    from esys.escript import printParallelThreadCounts
44    
45  #=========================================================  #=========================================================
46  #   some helpers:  #   some helpers:
47  #=========================================================  #=========================================================
48    def getEpsilon():
49         #     ------------------------------------------------------------------
50         #     Compute EPSILON, the machine precision.  The call to daxpp is
51         #     inTENded to fool compilers that use extra-length registers.
52         #     31 Map 1999: Hardwire EPSILON so the debugger can step thru easily.
53         #     ------------------------------------------------------------------
54         eps    = 2.**(-12)
55         p=2.
56         while p>1.:
57                eps/=2.
58                p=1.+eps
59         return eps*2.
60    
61    EPSILON=getEpsilon()
62    
63    def getTagNames(domain):
64        """
65        returns a list of the tag names used by the domain
66    
67        
68        @param domain: a domain object
69        @type domain: L{escript.Domain}
70        @return: a list of the tag name used by the domain.
71        @rtype: C{list} of C{str}
72        """
73        return [n.strip() for n in domain.showTagNames().split(",") ]
74    
75    def insertTagNames(domain,**kwargs):
76        """
77        inserts tag names into the domain
78    
79        @param domain: a domain object
80        @type domain: C{escript.Domain}
81        @keyword <tag name>: tag key assigned to <tag name>
82        @type <tag name>: C{int}
83        """
84        for  k in kwargs:
85             domain.setTagMap(k,kwargs[k])
86    
87    def insertTaggedValues(target,**kwargs):
88        """
89        inserts tagged values into the tagged using tag names
90    
91        @param target: data to be filled by tagged values
92        @type target: L{escript.Data}
93        @keyword <tag name>: value to be used for <tag name>
94        @type <tag name>: C{float} or {numarray.NumArray}
95        @return: C{target}
96        @rtype: L{escript.Data}
97        """
98        for k in kwargs:
99            target.setTaggedValue(k,kwargs[k])
100        return target
101    
102  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
103      """      """
104      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.
# Line 38  def saveVTK(filename,domain=None,**data) Line 107  def saveVTK(filename,domain=None,**data)
107    
108         tmp=Scalar(..)         tmp=Scalar(..)
109         v=Vector(..)         v=Vector(..)
110         saveVTK("solution.xml",temperature=tmp,velovity=v)         saveVTK("solution.xml",temperature=tmp,velocity=v)
111    
112      tmp and v are written into "solution.xml" where tmp is named "temperature" and v is named "velovity"      tmp and v are written into "solution.xml" where tmp is named "temperature" and v is named "velocity"
113    
114      @param filename: file name of the output file      @param filename: file name of the output file
115      @type filename: C{str}      @type filename: C{str}
# Line 50  def saveVTK(filename,domain=None,**data) Line 119  def saveVTK(filename,domain=None,**data)
119      @type <name>: L{Data} object.      @type <name>: L{Data} object.
120      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
121      """      """
122      if domain==None:      new_data={}
123         for i in data.keys():      for n,d in data.items():
124            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
125                fs=d.getFunctionSpace()
126                domain2=fs.getDomain()
127                if fs == escript.Solution(domain2):
128                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
129                elif fs == escript.ReducedSolution(domain2):
130                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
131                else:
132                   new_data[n]=d
133                if domain==None: domain=domain2
134      if domain==None:      if domain==None:
135          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
136      else:      domain.saveVTK(filename,new_data)
         domain.saveVTK(filename,data)  
137    
138  def saveDX(filename,domain=None,**data):  def saveDX(filename,domain=None,**data):
139      """      """
# Line 66  def saveDX(filename,domain=None,**data): Line 143  def saveDX(filename,domain=None,**data):
143    
144         tmp=Scalar(..)         tmp=Scalar(..)
145         v=Vector(..)         v=Vector(..)
146         saveDX("solution.dx",temperature=tmp,velovity=v)         saveDX("solution.dx",temperature=tmp,velocity=v)
147    
148      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity".      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velocity".
149    
150      @param filename: file name of the output file      @param filename: file name of the output file
151      @type filename: C{str}      @type filename: C{str}
# Line 78  def saveDX(filename,domain=None,**data): Line 155  def saveDX(filename,domain=None,**data):
155      @type <name>: L{Data} object.      @type <name>: L{Data} object.
156      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
157      """      """
158      if domain==None:      new_data={}
159         for i in data.keys():      for n,d in data.items():
160            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
161                fs=d.getFunctionSpace()
162                domain2=fs.getDomain()
163                if fs == escript.Solution(domain2):
164                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
165                elif fs == escript.ReducedSolution(domain2):
166                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
167                elif fs == escript.ContinuousFunction(domain2):
168                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
169                else:
170                   new_data[n]=d
171                if domain==None: domain=domain2
172      if domain==None:      if domain==None:
173          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
174      else:      domain.saveDX(filename,new_data)
         domain.saveDX(filename,data)  
175    
176  def kronecker(d=3):  def kronecker(d=3):
177     """     """
# Line 238  def inf(arg): Line 325  def inf(arg):
325  #=========================================================================  #=========================================================================
326  #   some little helpers  #   some little helpers
327  #=========================================================================  #=========================================================================
328  def pokeShape(arg):  def getRank(arg):
329        """
330        identifies the rank of its argument
331    
332        @param arg: a given object
333        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
334        @return: the rank of the argument
335        @rtype: C{int}
336        @raise TypeError: if type of arg cannot be processed
337        """
338    
339        if isinstance(arg,numarray.NumArray):
340            return arg.rank
341        elif isinstance(arg,escript.Data):
342            return arg.getRank()
343        elif isinstance(arg,float):
344            return 0
345        elif isinstance(arg,int):
346            return 0
347        elif isinstance(arg,Symbol):
348            return arg.getRank()
349        else:
350          raise TypeError,"getShape: cannot identify shape"
351    def getShape(arg):
352      """      """
353      identifies the shape of its argument      identifies the shape of its argument
354    
# Line 260  def pokeShape(arg): Line 370  def pokeShape(arg):
370      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
371          return arg.getShape()          return arg.getShape()
372      else:      else:
373        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
374    
375  def pokeDim(arg):  def pokeDim(arg):
376      """      """
# Line 283  def commonShape(arg0,arg1): Line 393  def commonShape(arg0,arg1):
393      """      """
394      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.
395    
396      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
397      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
398      @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.
399      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
400      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
401      """      """
402      sh0=pokeShape(arg0)      sh0=getShape(arg0)
403      sh1=pokeShape(arg1)      sh1=getShape(arg1)
404      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
405         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
406               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 440  def matchShape(arg0,arg1): Line 550  def matchShape(arg0,arg1):
550      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
551      @param arg1: a given object      @param arg1: a given object
552      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
553      @return: arg0 and arg1 where copies are returned when the shape has to be changed.      @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
554      @rtype: C{tuple}      @rtype: C{tuple}
555      """      """
556      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
557      sh0=pokeShape(arg0)      sh0=getShape(arg0)
558      sh1=pokeShape(arg1)      sh1=getShape(arg1)
559      if len(sh0)<len(sh):      if len(sh0)<len(sh):
560         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
561      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
# Line 573  class Symbol(object): Line 683  class Symbol(object):
683            if isinstance(a,Symbol):            if isinstance(a,Symbol):
684               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
685            else:            else:
686                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
687                if len(s)>0:                if len(s)>0:
688                   out.append(numarray.zeros(s),numarray.Float64)                   out.append(numarray.zeros(s),numarray.Float64)
689                else:                else:
# Line 1310  def whereNonZero(arg,tol=0.): Line 1420  def whereNonZero(arg,tol=0.):
1420     else:     else:
1421        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1422    
1423    def erf(arg):
1424       """
1425       returns erf of argument arg
1426    
1427       @param arg: argument
1428       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1429       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1430       @raises TypeError: if the type of the argument is not expected.
1431       """
1432       if isinstance(arg,escript.Data):
1433          return arg._erf()
1434       else:
1435          raise TypeError,"erf: Unknown argument type."
1436    
1437  def sin(arg):  def sin(arg):
1438     """     """
1439     returns sine of argument arg     returns sine of argument arg
# Line 2930  def trace(arg,axis_offset=0): Line 3054  def trace(arg,axis_offset=0):
3054    
3055     @param arg: argument     @param arg: argument
3056     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3057     @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
3058                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3059     @type axis_offset: C{int}     @type axis_offset: C{int}
3060     @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.
3061     @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 2963  def trace(arg,axis_offset=0): Line 3087  def trace(arg,axis_offset=0):
3087        s=list(arg.getShape())                s=list(arg.getShape())        
3088        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3089          raise ValueError,"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)
3090        return arg._matrixtrace(axis_offset)        return arg._trace(axis_offset)
3091     elif isinstance(arg,float):     elif isinstance(arg,float):
3092        raise TypeError,"illegal argument type float."        raise TypeError,"illegal argument type float."
3093     elif isinstance(arg,int):     elif isinstance(arg,int):
# Line 2982  class Trace_Symbol(DependendSymbol): Line 3106  class Trace_Symbol(DependendSymbol):
3106        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3107        @param arg: argument of function        @param arg: argument of function
3108        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3109        @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
3110                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3111        @type axis_offset: C{int}        @type axis_offset: C{int}
3112        """        """
3113        if arg.getRank()<2:        if arg.getRank()<2:
# Line 3049  class Trace_Symbol(DependendSymbol): Line 3173  class Trace_Symbol(DependendSymbol):
3173    
3174  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3175     """     """
3176     returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.     returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3177    
3178     @param arg: argument     @param arg: argument
3179     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3180     @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.     @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3181                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3182     @type axis_offset: C{int}     @type axis_offset: C{int}
3183     @return: transpose of arg     @return: transpose of arg
3184     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
# Line 3092  class Transpose_Symbol(DependendSymbol): Line 3216  class Transpose_Symbol(DependendSymbol):
3216    
3217        @param arg: argument of function        @param arg: argument of function
3218        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3219        @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.        @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3220                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3221        @type axis_offset: C{int}        @type axis_offset: C{int}
3222        """        """
3223        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3224        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3225          raise ValueError,"axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3226        s=arg.getShape()        s=arg.getShape()
3227        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3228    
# Line 3153  class Transpose_Symbol(DependendSymbol): Line 3277  class Transpose_Symbol(DependendSymbol):
3277           return identity(self.getShape())           return identity(self.getShape())
3278        else:        else:
3279           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3280    
3281    def swap_axes(arg,axis0=0,axis1=1):
3282       """
3283       returns the swap of arg by swaping the components axis0 and axis1
3284    
3285       @param arg: argument
3286       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3287       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3288       @type axis0: C{int}
3289       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3290       @type axis1: C{int}
3291       @return: C{arg} with swaped components
3292       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3293       """
3294       if axis0 > axis1:
3295          axis0,axis1=axis1,axis0
3296       if isinstance(arg,numarray.NumArray):
3297          return numarray.swapaxes(arg,axis0,axis1)
3298       elif isinstance(arg,escript.Data):
3299          return arg._swap_axes(axis0,axis1)
3300       elif isinstance(arg,float):
3301          raise TyepError,"float argument is not supported."
3302       elif isinstance(arg,int):
3303          raise TyepError,"int argument is not supported."
3304       elif isinstance(arg,Symbol):
3305          return SwapAxes_Symbol(arg,axis0,axis1)
3306       else:
3307          raise TypeError,"Unknown argument type."
3308    
3309    class SwapAxes_Symbol(DependendSymbol):
3310       """
3311       L{Symbol} representing the result of the swap function
3312       """
3313       def __init__(self,arg,axis0=0,axis1=1):
3314          """
3315          initialization of swap L{Symbol} with argument arg
3316    
3317          @param arg: argument
3318          @type arg: L{Symbol}.
3319          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3320          @type axis0: C{int}
3321          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3322          @type axis1: C{int}
3323          """
3324          if arg.getRank()<2:
3325             raise ValueError,"argument must have at least rank 2."
3326          if axis0<0 or axis0>arg.getRank()-1:
3327             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3328          if axis1<0 or axis1>arg.getRank()-1:
3329             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3330          if axis0 == axis1:
3331             raise ValueError,"axis indices must be different."
3332          if axis0 > axis1:
3333             axis0,axis1=axis1,axis0
3334          s=arg.getShape()
3335          s_out=[]
3336          for i in range(len(s)):
3337             if i == axis0:
3338                s_out.append(s[axis1])
3339             elif i == axis1:
3340                s_out.append(s[axis0])
3341             else:
3342                s_out.append(s[i])
3343          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3344    
3345       def getMyCode(self,argstrs,format="escript"):
3346          """
3347          returns a program code that can be used to evaluate the symbol.
3348    
3349          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3350          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3351          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3352          @type format: C{str}
3353          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3354          @rtype: C{str}
3355          @raise NotImplementedError: if the requested format is not available
3356          """
3357          if format=="escript" or format=="str"  or format=="text":
3358             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3359          else:
3360             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3361    
3362       def substitute(self,argvals):
3363          """
3364          assigns new values to symbols in the definition of the symbol.
3365          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3366    
3367          @param argvals: new values assigned to symbols
3368          @type argvals: C{dict} with keywords of type L{Symbol}.
3369          @return: result of the substitution process. Operations are executed as much as possible.
3370          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3371          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3372          """
3373          if argvals.has_key(self):
3374             arg=argvals[self]
3375             if self.isAppropriateValue(arg):
3376                return arg
3377             else:
3378                raise TypeError,"%s: new value is not appropriate."%str(self)
3379          else:
3380             arg=self.getSubstitutedArguments(argvals)
3381             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3382    
3383       def diff(self,arg):
3384          """
3385          differential of this object
3386    
3387          @param arg: the derivative is calculated with respect to arg
3388          @type arg: L{escript.Symbol}
3389          @return: derivative with respect to C{arg}
3390          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3391          """
3392          if arg==self:
3393             return identity(self.getShape())
3394          else:
3395             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3396    
3397  def symmetric(arg):  def symmetric(arg):
3398      """      """
3399      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
# Line 3528  class Add_Symbol(DependendSymbol): Line 3769  class Add_Symbol(DependendSymbol):
3769         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3770         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3771         """         """
3772         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3773         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3774         if not sh0==sh1:         if not sh0==sh1:
3775            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3776         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 3572  class Add_Symbol(DependendSymbol): Line 3813  class Add_Symbol(DependendSymbol):
3813              raise TypeError,"%s: new value is not appropriate."%str(self)              raise TypeError,"%s: new value is not appropriate."%str(self)
3814        else:        else:
3815           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
3816           return add(args[0],args[1])           out=add(args[0],args[1])
3817             return out
3818    
3819     def diff(self,arg):     def diff(self,arg):
3820        """        """
# Line 3603  def mult(arg0,arg1): Line 3845  def mult(arg0,arg1):
3845         """         """
3846         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3847         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3848            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3849         else:         else:
3850            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3851                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3627  class Mult_Symbol(DependendSymbol): Line 3869  class Mult_Symbol(DependendSymbol):
3869         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3870         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3871         """         """
3872         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3873         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3874         if not sh0==sh1:         if not sh0==sh1:
3875            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3876         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 3703  def quotient(arg0,arg1): Line 3945  def quotient(arg0,arg1):
3945         """         """
3946         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3947         if testForZero(args[0]):         if testForZero(args[0]):
3948            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3949         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3950            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3951               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3732  class Quotient_Symbol(DependendSymbol): Line 3974  class Quotient_Symbol(DependendSymbol):
3974         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3975         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3976         """         """
3977         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3978         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3979         if not sh0==sh1:         if not sh0==sh1:
3980            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3981         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 3809  def power(arg0,arg1): Line 4051  def power(arg0,arg1):
4051         """         """
4052         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4053         if testForZero(args[0]):         if testForZero(args[0]):
4054            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4055         elif testForZero(args[1]):         elif testForZero(args[1]):
4056            return numarray.ones(pokeShape(args[1]),numarray.Float64)            return numarray.ones(getShape(args[1]),numarray.Float64)
4057         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4058            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4059         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 3834  class Power_Symbol(DependendSymbol): Line 4076  class Power_Symbol(DependendSymbol):
4076         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4077         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4078         """         """
4079         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4080         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4081         if not sh0==sh1:         if not sh0==sh1:
4082            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4083         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3934  def minimum(*args): Line 4176  def minimum(*args):
4176            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4177      return out      return out
4178    
4179  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4180      """      """
4181      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4182    
4183      @param arg: argument      @param arg: argument
4184      @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}
4185      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4186      @type minval: C{float}      @type minval: C{float} or C{None}
4187      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4188      @type maxval: C{float}      @type maxval: C{float} or C{None}
4189      @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
4190               less then maxval are unchanged.               less then maxval are unchanged.
4191      @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
4192      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4193      """      """
4194      if minval>maxval:      if not minval==None and not maxval==None:
4195         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4196      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4197        if minval == None:
4198            tmp=arg
4199        else:
4200            tmp=maximum(minval,arg)
4201        if maxval == None:
4202            return tmp
4203        else:
4204            return minimum(tmp,maxval)
4205    
4206        
4207  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3972  def inner(arg0,arg1): Line 4222  def inner(arg0,arg1):
4222      @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
4223      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4224      """      """
4225      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4226      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4227      if not sh0==sh1:      if not sh0==sh1:
4228          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4229      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
# Line 4024  def matrix_mult(arg0,arg1): Line 4274  def matrix_mult(arg0,arg1):
4274      @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
4275      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4276      """      """
4277      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4278      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4279      if not len(sh0)==2 :      if not len(sh0)==2 :
4280          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4281      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
# Line 4073  def tensor_mult(arg0,arg1): Line 4323  def tensor_mult(arg0,arg1):
4323      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4324      @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
4325      """      """
4326      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4327      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4328      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4329         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4330      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):
# Line 4148  class GeneralTensorProduct_Symbol(Depend Line 4398  class GeneralTensorProduct_Symbol(Depend
4398         @raise ValueError: illegal dimension         @raise ValueError: illegal dimension
4399         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4400         """         """
4401         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4402         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4403         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4404         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4405         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 4196  class GeneralTensorProduct_Symbol(Depend Line 4446  class GeneralTensorProduct_Symbol(Depend
4446           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4447           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4448    
4449  def escript_generalTensorProductNew(arg0,arg1,axis_offset):  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4450      "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!!!"
4451      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
     shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
     shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  
     shape10=arg1.getShape()[:axis_offset]  
     shape1=arg1.getShape()[axis_offset:]  
     if not shape01==shape10:  
         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)  
     # Figure out which functionspace to use (look at where operator+ is defined maybe in BinaryOp.h to get the logic for this)  
     # fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
     out=GeneralTensorProduct(arg0, arg1, axis_offset)  
     return out  
   
 def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  
     "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"  
     # calculate the return shape:  
     shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
     shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  
     shape10=arg1.getShape()[:axis_offset]  
     shape1=arg1.getShape()[axis_offset:]  
     if not shape01==shape10:  
         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)  
   
     # 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  
   
4452    
4453  def transposed_matrix_mult(arg0,arg1):  def transposed_matrix_mult(arg0,arg1):
4454      """      """
# Line 4275  def transposed_matrix_mult(arg0,arg1): Line 4472  def transposed_matrix_mult(arg0,arg1):
4472      @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
4473      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4474      """      """
4475      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4476      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4477      if not len(sh0)==2 :      if not len(sh0)==2 :
4478          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4479      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
# Line 4320  def transposed_tensor_mult(arg0,arg1): Line 4517  def transposed_tensor_mult(arg0,arg1):
4517      @return: the tensor product of tarnsposed of arg0 and arg1 at each data point      @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4518      @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
4519      """      """
4520      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4521      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4522      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4523         return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)         return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4524      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):
# Line 4398  class GeneralTransposedTensorProduct_Sym Line 4595  class GeneralTransposedTensorProduct_Sym
4595         @raise ValueError: inconsistent dimensions of arguments.         @raise ValueError: inconsistent dimensions of arguments.
4596         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4597         """         """
4598         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4599         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4600         sh01=sh_arg0[:axis_offset]         sh01=sh_arg0[:axis_offset]
4601         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
4602         sh0=sh_arg0[axis_offset:]         sh0=sh_arg0[axis_offset:]
# Line 4448  class GeneralTransposedTensorProduct_Sym Line 4645  class GeneralTransposedTensorProduct_Sym
4645    
4646  def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct  def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4647      "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!!!"
4648      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
     shape01=arg0.getShape()[:axis_offset]  
     shape10=arg1.getShape()[:axis_offset]  
     shape0=arg0.getShape()[axis_offset:]  
     shape1=arg1.getShape()[axis_offset:]  
     if not shape01==shape10:  
         raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)  
   
     # 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(i01+i0))*arg1.__getitem__(tuple(i01+i1))  
          out.__setitem__(tuple(i0+i1),s)  
     return out  
   
4649    
4650  def matrix_transposed_mult(arg0,arg1):  def matrix_transposed_mult(arg0,arg1):
4651      """      """
# Line 4507  def matrix_transposed_mult(arg0,arg1): Line 4665  def matrix_transposed_mult(arg0,arg1):
4665      @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
4666      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4667      """      """
4668      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4669      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4670      if not len(sh0)==2 :      if not len(sh0)==2 :
4671          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4672      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
# Line 4543  def tensor_transposed_mult(arg0,arg1): Line 4701  def tensor_transposed_mult(arg0,arg1):
4701      @return: the tensor product of tarnsposed of arg0 and arg1 at each data point      @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4702      @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
4703      """      """
4704      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4705      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4706      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4707         return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)         return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4708      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):
# Line 4621  class GeneralTensorTransposedProduct_Sym Line 4779  class GeneralTensorTransposedProduct_Sym
4779         @raise ValueError: inconsistent dimensions of arguments.         @raise ValueError: inconsistent dimensions of arguments.
4780         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4781         """         """
4782         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4783         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4784         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4785         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4786         sh10=sh_arg1[len(sh_arg1)-axis_offset:]         sh10=sh_arg1[len(sh_arg1)-axis_offset:]
# Line 4671  class GeneralTensorTransposedProduct_Sym Line 4829  class GeneralTensorTransposedProduct_Sym
4829    
4830  def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct  def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4831      "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!!!"
4832      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
     shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  
     shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
     shape10=arg1.getShape()[arg1.getRank()-axis_offset:]  
     shape1=arg1.getShape()[:arg1.getRank()-axis_offset]  
     if not shape01==shape10:  
         raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)  
   
     # 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(i1+i01))  
          out.__setitem__(tuple(i0+i1),s)  
     return out  
   
4833    
4834  #=========================================================  #=========================================================
4835  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
# Line 4836  def integrate(arg,where=None): Line 4955  def integrate(arg,where=None):
4955         else:         else:
4956            return arg._integrate()            return arg._integrate()
4957      else:      else:
4958        raise TypeError,"integrate: Unknown argument type."         arg2=escript.Data(arg,where)
4959           if arg2.getRank()==0:
4960              return arg2._integrate()[0]
4961           else:
4962              return arg2._integrate()
4963    
4964  class Integrate_Symbol(DependendSymbol):  class Integrate_Symbol(DependendSymbol):
4965     """     """
# Line 4908  class Integrate_Symbol(DependendSymbol): Line 5031  class Integrate_Symbol(DependendSymbol):
5031    
5032  def interpolate(arg,where):  def interpolate(arg,where):
5033      """      """
5034      interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where. If the argument C{arg} has the requested function space
5035        C{where} no interpolation is performed and C{arg} is returned.
5036    
5037      @param arg: interpolant      @param arg: interpolant
5038      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
# Line 4919  def interpolate(arg,where): Line 5043  def interpolate(arg,where):
5043      """      """
5044      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5045         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
5046        elif isinstance(arg,escript.Data):
5047           if where == arg.getFunctionSpace():
5048              return arg
5049           else:
5050              return escript.Data(arg,where)
5051      else:      else:
5052         return escript.Data(arg,where)         return escript.Data(arg,where)
5053    
# Line 5037  def L2(arg): Line 5166  def L2(arg):
5166      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5167      """      """
5168      return sqrt(integrate(inner(arg,arg)))      return sqrt(integrate(inner(arg,arg)))
5169    
5170    def getClosestValue(arg,origin=0):
5171        """
5172        returns the value in arg which is closest to origin
5173        
5174        @param arg: function
5175        @type arg: L{escript.Data} or L{Symbol}
5176        @param origin: reference value
5177        @type origin: C{float} or L{escript.Data}
5178        @return: value in arg closest to origin
5179        @rtype: L{numarray.NumArray} or L{Symbol}
5180        """
5181        return arg.getValueOfGlobalDataPoint(*(length(arg-origin).minGlobalDataPoint()))
5182    
5183    def normalize(arg,zerolength=0):
5184        """
5185        returns normalized version of arg (=arg/length(arg))
5186        
5187        @param arg: function
5188        @type arg: L{escript.Data} or L{Symbol}
5189        @param zerolength: realitive tolerance for arg == 0.
5190        @type zerolength: C{float}
5191        @return: normalized arg where arg is non zero and zero elsewhere
5192        @rtype: L{escript.Data} or L{Symbol}
5193        """
5194        l=length(arg)
5195        m=whereZero(l,zerolength*Lsup(l))
5196        mm=1-m
5197        return arg*(mm/(l*mm+m))
5198    
5199  #=============================  #=============================
5200  #  #
5201    
# Line 5046  def reorderComponents(arg,index): Line 5205  def reorderComponents(arg,index):
5205    
5206      """      """
5207      raise NotImplementedError      raise NotImplementedError
 #  
 # $Log: util.py,v $  
 # Revision 1.14.2.16  2005/10/19 06:09:57  gross  
 # new versions of saveVTK and saveDX which allow to write several Data objects into a single file  
 #  
 # Revision 1.14.2.15  2005/10/19 03:25:27  jgs  
 # numarray uses log10 for log, and log for ln - log()/ln() updated to reflect this  
 #  
 # Revision 1.14.2.14  2005/10/19 02:10:23  jgs  
 # fixed ln() to call Data::ln  
 #  
 # Revision 1.14.2.13  2005/09/12 03:32:14  gross  
 # test_visualiztion has been aded to mk  
 #  
 # Revision 1.14.2.12  2005/09/09 01:56:24  jgs  
 # added implementations of acos asin atan sinh cosh tanh asinh acosh atanh  
 # and some associated testing  
 #  
 # Revision 1.14.2.11  2005/09/08 08:28:39  gross  
 # some cleanup in savevtk  
 #  
 # Revision 1.14.2.10  2005/09/08 00:25:32  gross  
 # test for finley mesh generators added  
 #  
 # Revision 1.14.2.9  2005/09/07 10:32:05  gross  
 # Symbols removed from util and put into symmbols.py.  
 #  
 # Revision 1.14.2.8  2005/08/26 05:06:37  cochrane  
 # Corrected errors in docstrings.  Improved output formatting of docstrings.  
 # Other minor improvements to code and docs (eg spelling etc).  
 #  
 # Revision 1.14.2.7  2005/08/26 04:45:40  cochrane  
 # Fixed and tidied markup and docstrings.  Some *Symbol classes were defined  
 # as functions, so changed them to classes (hopefully this was the right thing  
 # to do).  
 #  
 # Revision 1.14.2.6  2005/08/26 04:30:13  gross  
 # gneric unit testing for linearPDE  
 #  
 # Revision 1.14.2.5  2005/08/24 02:02:52  gross  
 # jump function added  
 #  
 # Revision 1.14.2.4  2005/08/18 04:39:32  gross  
 # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now  
 #  
 # Revision 1.14.2.3  2005/08/03 09:55:33  gross  
 # ContactTest is passing now./mk install!  
 #  
 # Revision 1.14.2.2  2005/08/02 03:15:14  gross  
 # bug inb trace fixed!  
 #  
 # Revision 1.14.2.1  2005/07/29 07:10:28  gross  
 # new functions in util and a new pde type in linearPDEs  
 #  
 # Revision 1.2.2.21  2005/07/28 04:19:23  gross  
 # new functions maximum and minimum introduced.  
 #  
 # Revision 1.2.2.20  2005/07/25 01:26:27  gross  
 # bug in inner fixed  
 #  
 # Revision 1.2.2.19  2005/07/21 04:01:28  jgs  
 # minor comment fixes  
 #  
 # Revision 1.2.2.18  2005/07/21 01:02:43  jgs  
 # commit ln() updates to development branch version  
 #  
 # Revision 1.12  2005/07/20 06:14:58  jgs  
 # added ln(data) style wrapper for data.ln() - also added corresponding  
 # implementation of Ln_Symbol class (not sure if this is right though)  
 #  
 # Revision 1.11  2005/07/08 04:07:35  jgs  
 # Merge of development branch back to main trunk on 2005-07-08  
 #  
 # Revision 1.10  2005/06/09 05:37:59  jgs  
 # Merge of development branch back to main trunk on 2005-06-09  
 #  
 # Revision 1.2.2.17  2005/07/07 07:28:58  gross  
 # some stuff added to util.py to improve functionality  
 #  
 # Revision 1.2.2.16  2005/06/30 01:53:55  gross  
 # a bug in coloring fixed  
 #  
 # Revision 1.2.2.15  2005/06/29 02:36:43  gross  
 # Symbols have been introduced and some function clarified. needs much more work  
 #  
 # Revision 1.2.2.14  2005/05/20 04:05:23  gross  
 # some work on a darcy flow started  
 #  
 # Revision 1.2.2.13  2005/03/16 05:17:58  matt  
 # Implemented unit(idx, dim) to create cartesian unit basis vectors to  
 # complement kronecker(dim) function.  
 #  
 # Revision 1.2.2.12  2005/03/10 08:14:37  matt  
 # Added non-member Linf utility function to complement Data::Linf().  
 #  
 # Revision 1.2.2.11  2005/02/17 05:53:25  gross  
 # some bug in saveDX fixed: in fact the bug was in  
 # DataC/getDataPointShape  
 #  
 # Revision 1.2.2.10  2005/01/11 04:59:36  gross  
 # automatic interpolation in integrate switched off  
 #  
 # Revision 1.2.2.9  2005/01/11 03:38:13  gross  
 # Bug in Data.integrate() fixed for the case of rank 0. The problem is not totallly resolved as the method should return a scalar rather than a numarray object in the case of rank 0. This problem is fixed by the util.integrate wrapper.  
 #  
 # Revision 1.2.2.8  2005/01/05 04:21:41  gross  
 # FunctionSpace checking/matchig in slicing added  
 #  
 # Revision 1.2.2.7  2004/12/29 05:29:59  gross  
 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()  
 #  
 # Revision 1.2.2.6  2004/12/24 06:05:41  gross  
 # some changes in linearPDEs to add AdevectivePDE  
 #  
 # Revision 1.2.2.5  2004/12/17 00:06:53  gross  
 # mk sets ESYS_ROOT is undefined  
 #  
 # Revision 1.2.2.4  2004/12/07 03:19:51  gross  
 # options for GMRES and PRES20 added  
 #  
 # Revision 1.2.2.3  2004/12/06 04:55:18  gross  
 # function wraper extended  
 #  
 # Revision 1.2.2.2  2004/11/22 05:44:07  gross  
 # a few more unitary functions have been added but not implemented in Data yet  
 #  
 # Revision 1.2.2.1  2004/11/12 06:58:15  gross  
 # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry  
 #  
 # Revision 1.2  2004/10/27 00:23:36  jgs  
 # fixed minor syntax error  
 #  
 # Revision 1.1.1.1  2004/10/26 06:53:56  jgs  
 # initial import of project esys2  
 #  
 # Revision 1.1.2.3  2004/10/26 06:43:48  jgs  
 # committing Lutz's and Paul's changes to brach jgs  
 #  
 # Revision 1.1.4.1  2004/10/20 05:32:51  cochrane  
 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.  
 #  
 # Revision 1.1  2004/08/05 03:58:27  gross  
 # Bug in Assemble_NodeCoordinates fixed  
 #  
 #  
   
 # vim: expandtab shiftwidth=4:  

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

  ViewVC Help
Powered by ViewVC 1.1.26