/[escript]/trunk-mpi-branch/escript/py_src/util.py
ViewVC logotype

Diff of /trunk-mpi-branch/escript/py_src/util.py

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

trunk/escript/py_src/util.py revision 775 by ksteube, Mon Jul 10 04:00:08 2006 UTC trunk-mpi-branch/escript/py_src/util.py revision 1295 by ksteube, Mon Sep 10 06:07:09 2007 UTC
# Line 26  import math Line 26  import math
26  import numarray  import numarray
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:  #   some helpers:
34  #=========================================================  #=========================================================
35    def getTagNames(domain):
36        """
37        returns a list of the tag names used by the domain
38    
39        
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 insertTagNames(domain,**kwargs):
48        """
49        inserts tag names into the domain
50    
51        @param domain: a domain object
52        @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        """
61        inserts tagged values into the tagged using tag names
62    
63        @param target: data to be filled by tagged values
64        @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  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
75      """      """
76      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 50  def saveVTK(filename,domain=None,**data) Line 91  def saveVTK(filename,domain=None,**data)
91      @type <name>: L{Data} object.      @type <name>: L{Data} object.
92      @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.
93      """      """
94      if domain==None:      new_data={}
95         for i in data.keys():      for n,d in data.items():
96            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
97                fs=d.getFunctionSpace()
98                domain2=fs.getDomain()
99                if fs == escript.Solution(domain2):
100                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
101                elif fs == escript.ReducedSolution(domain2):
102                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
103                else:
104                   new_data[n]=d
105                if domain==None: domain=domain2
106      if domain==None:      if domain==None:
107          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
108      else:      domain.saveVTK(filename,new_data)
         domain.saveVTK(filename,data)  
109    
110  def saveDX(filename,domain=None,**data):  def saveDX(filename,domain=None,**data):
111      """      """
# Line 78  def saveDX(filename,domain=None,**data): Line 127  def saveDX(filename,domain=None,**data):
127      @type <name>: L{Data} object.      @type <name>: L{Data} object.
128      @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.
129      """      """
130      if domain==None:      new_data={}
131         for i in data.keys():      for n,d in data.items():
132            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
133                fs=d.getFunctionSpace()
134                domain2=fs.getDomain()
135                if fs == escript.Solution(domain2):
136                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
137                elif fs == escript.ReducedSolution(domain2):
138                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
139                elif fs == escript.ContinuousFunction(domain2):
140                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
141                else:
142                   new_data[n]=d
143                if domain==None: domain=domain2
144      if domain==None:      if domain==None:
145          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
146      else:      domain.saveDX(filename,new_data)
         domain.saveDX(filename,data)  
147    
148  def kronecker(d=3):  def kronecker(d=3):
149     """     """
# Line 238  def inf(arg): Line 297  def inf(arg):
297  #=========================================================================  #=========================================================================
298  #   some little helpers  #   some little helpers
299  #=========================================================================  #=========================================================================
300  def pokeShape(arg):  def getRank(arg):
301        """
302        identifies the rank of its argument
303    
304        @param arg: a given object
305        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
306        @return: the rank of the argument
307        @rtype: C{int}
308        @raise TypeError: if type of arg cannot be processed
309        """
310    
311        if isinstance(arg,numarray.NumArray):
312            return arg.rank
313        elif isinstance(arg,escript.Data):
314            return arg.getRank()
315        elif isinstance(arg,float):
316            return 0
317        elif isinstance(arg,int):
318            return 0
319        elif isinstance(arg,Symbol):
320            return arg.getRank()
321        else:
322          raise TypeError,"getShape: cannot identify shape"
323    def getShape(arg):
324      """      """
325      identifies the shape of its argument      identifies the shape of its argument
326    
# Line 260  def pokeShape(arg): Line 342  def pokeShape(arg):
342      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
343          return arg.getShape()          return arg.getShape()
344      else:      else:
345        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
346    
347  def pokeDim(arg):  def pokeDim(arg):
348      """      """
# Line 283  def commonShape(arg0,arg1): Line 365  def commonShape(arg0,arg1):
365      """      """
366      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.
367    
368      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
369      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
370      @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.
371      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
372      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
373      """      """
374      sh0=pokeShape(arg0)      sh0=getShape(arg0)
375      sh1=pokeShape(arg1)      sh1=getShape(arg1)
376      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
377         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
378               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 522  def matchShape(arg0,arg1):
522      @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}
523      @param arg1: a given object      @param arg1: a given object
524      @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}
525      @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.
526      @rtype: C{tuple}      @rtype: C{tuple}
527      """      """
528      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
529      sh0=pokeShape(arg0)      sh0=getShape(arg0)
530      sh1=pokeShape(arg1)      sh1=getShape(arg1)
531      if len(sh0)<len(sh):      if len(sh0)<len(sh):
532         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
533      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
# Line 573  class Symbol(object): Line 655  class Symbol(object):
655            if isinstance(a,Symbol):            if isinstance(a,Symbol):
656               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
657            else:            else:
658                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
659                if len(s)>0:                if len(s)>0:
660                   out.append(numarray.zeros(s),numarray.Float64)                   out.append(numarray.zeros(s),numarray.Float64)
661                else:                else:
# Line 1310  def whereNonZero(arg,tol=0.): Line 1392  def whereNonZero(arg,tol=0.):
1392     else:     else:
1393        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1394    
1395    def erf(arg):
1396       """
1397       returns erf of argument arg
1398    
1399       @param arg: argument
1400       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1401       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1402       @raises TypeError: if the type of the argument is not expected.
1403       """
1404       if isinstance(arg,escript.Data):
1405          return arg._erf()
1406       else:
1407          raise TypeError,"erf: Unknown argument type."
1408    
1409  def sin(arg):  def sin(arg):
1410     """     """
1411     returns sine of argument arg     returns sine of argument arg
# Line 2930  def trace(arg,axis_offset=0): Line 3026  def trace(arg,axis_offset=0):
3026    
3027     @param arg: argument     @param arg: argument
3028     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3029     @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
3030                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3031     @type axis_offset: C{int}     @type axis_offset: C{int}
3032     @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.
3033     @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 2939  def trace(arg,axis_offset=0): Line 3035  def trace(arg,axis_offset=0):
3035     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
3036        sh=arg.shape        sh=arg.shape
3037        if len(sh)<2:        if len(sh)<2:
3038          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3039        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
3040          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
3041        s1=1        s1=1
3042        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
3043        s2=1        s2=1
3044        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3045        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
3046          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)
3047        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))
3048        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
3049        for i1 in range(s1):        for i1 in range(s1):
# Line 2956  def trace(arg,axis_offset=0): Line 3052  def trace(arg,axis_offset=0):
3052        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3053        return out        return out
3054     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3055        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
3056            raise ValueError,"rank of argument must be greater than 1"
3057          if axis_offset<0 or axis_offset>arg.getRank()-2:
3058            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3059          s=list(arg.getShape())        
3060          if not s[axis_offset] == s[axis_offset+1]:
3061            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3062          return arg._trace(axis_offset)
3063     elif isinstance(arg,float):     elif isinstance(arg,float):
3064        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
3065     elif isinstance(arg,int):     elif isinstance(arg,int):
3066        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
3067     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3068        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
3069     else:     else:
3070        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
3071    
 def escript_trace(arg,axis_offset):  
       "arg is a Data object"  
       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=arg._matrixtrace(axis_offset)  
       return out  
3072  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
3073     """     """
3074     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 2986  class Trace_Symbol(DependendSymbol): Line 3078  class Trace_Symbol(DependendSymbol):
3078        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3079        @param arg: argument of function        @param arg: argument of function
3080        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3081        @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
3082                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3083        @type axis_offset: C{int}        @type axis_offset: C{int}
3084        """        """
3085        if arg.getRank()<2:        if arg.getRank()<2:
3086          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3087        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3088          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
3089        s=list(arg.getShape())                s=list(arg.getShape())        
3090        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3091          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)
3092        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())
3093    
3094     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 3053  class Trace_Symbol(DependendSymbol): Line 3145  class Trace_Symbol(DependendSymbol):
3145    
3146  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3147     """     """
3148     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.
3149    
3150     @param arg: argument     @param arg: argument
3151     @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}
3152     @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.
3153                         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.
3154     @type axis_offset: C{int}     @type axis_offset: C{int}
3155     @return: transpose of arg     @return: transpose of arg
3156     @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 3067  def transpose(arg,axis_offset=None): Line 3159  def transpose(arg,axis_offset=None):
3159        if axis_offset==None: axis_offset=int(arg.rank/2)        if axis_offset==None: axis_offset=int(arg.rank/2)
3160        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3161     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3162        if axis_offset==None: axis_offset=int(arg.getRank()/2)        r=arg.getRank()
3163        return escript_transpose(arg,axis_offset)        if axis_offset==None: axis_offset=int(r/2)
3164          if axis_offset<0 or axis_offset>r:
3165            raise ValueError,"axis_offset must be between 0 and %s"%r
3166          return arg._transpose(axis_offset)
3167     elif isinstance(arg,float):     elif isinstance(arg,float):
3168        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3169          raise ValueError,"transpose: axis_offset must be 0 for float argument"          raise ValueError,"axis_offset must be 0 for float argument"
3170        return arg        return arg
3171     elif isinstance(arg,int):     elif isinstance(arg,int):
3172        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3173          raise ValueError,"transpose: axis_offset must be 0 for int argument"          raise ValueError,"axis_offset must be 0 for int argument"
3174        return float(arg)        return float(arg)
3175     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3176        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3177        return Transpose_Symbol(arg,axis_offset)        return Transpose_Symbol(arg,axis_offset)
3178     else:     else:
3179        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"Unknown argument type."
3180    
 def escript_transpose(arg,axis_offset):  
       "arg is a Data object"  
       r=arg.getRank()  
       if axis_offset<0 or axis_offset>r:  
         raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r  
       out=arg._transpose(axis_offset)  
       return out  
3181  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3182     """     """
3183     L{Symbol} representing the result of the transpose function     L{Symbol} representing the result of the transpose function
# Line 3100  class Transpose_Symbol(DependendSymbol): Line 3188  class Transpose_Symbol(DependendSymbol):
3188    
3189        @param arg: argument of function        @param arg: argument of function
3190        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3191        @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.
3192                         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.
3193        @type axis_offset: C{int}        @type axis_offset: C{int}
3194        """        """
3195        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3196        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3197          raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3198        s=arg.getShape()        s=arg.getShape()
3199        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())
3200    
# Line 3161  class Transpose_Symbol(DependendSymbol): Line 3249  class Transpose_Symbol(DependendSymbol):
3249           return identity(self.getShape())           return identity(self.getShape())
3250        else:        else:
3251           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3252    
3253    def swap_axes(arg,axis0=0,axis1=1):
3254       """
3255       returns the swap of arg by swaping the components axis0 and axis1
3256    
3257       @param arg: argument
3258       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3259       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3260       @type axis0: C{int}
3261       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3262       @type axis1: C{int}
3263       @return: C{arg} with swaped components
3264       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3265       """
3266       if axis0 > axis1:
3267          axis0,axis1=axis1,axis0
3268       if isinstance(arg,numarray.NumArray):
3269          return numarray.swapaxes(arg,axis0,axis1)
3270       elif isinstance(arg,escript.Data):
3271          return arg._swap_axes(axis0,axis1)
3272       elif isinstance(arg,float):
3273          raise TyepError,"float argument is not supported."
3274       elif isinstance(arg,int):
3275          raise TyepError,"int argument is not supported."
3276       elif isinstance(arg,Symbol):
3277          return SwapAxes_Symbol(arg,axis0,axis1)
3278       else:
3279          raise TypeError,"Unknown argument type."
3280    
3281    class SwapAxes_Symbol(DependendSymbol):
3282       """
3283       L{Symbol} representing the result of the swap function
3284       """
3285       def __init__(self,arg,axis0=0,axis1=1):
3286          """
3287          initialization of swap L{Symbol} with argument arg
3288    
3289          @param arg: argument
3290          @type arg: L{Symbol}.
3291          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3292          @type axis0: C{int}
3293          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3294          @type axis1: C{int}
3295          """
3296          if arg.getRank()<2:
3297             raise ValueError,"argument must have at least rank 2."
3298          if axis0<0 or axis0>arg.getRank()-1:
3299             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3300          if axis1<0 or axis1>arg.getRank()-1:
3301             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3302          if axis0 == axis1:
3303             raise ValueError,"axis indices must be different."
3304          if axis0 > axis1:
3305             axis0,axis1=axis1,axis0
3306          s=arg.getShape()
3307          s_out=[]
3308          for i in range(len(s)):
3309             if i == axis0:
3310                s_out.append(s[axis1])
3311             elif i == axis1:
3312                s_out.append(s[axis0])
3313             else:
3314                s_out.append(s[i])
3315          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3316    
3317       def getMyCode(self,argstrs,format="escript"):
3318          """
3319          returns a program code that can be used to evaluate the symbol.
3320    
3321          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3322          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3323          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3324          @type format: C{str}
3325          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3326          @rtype: C{str}
3327          @raise NotImplementedError: if the requested format is not available
3328          """
3329          if format=="escript" or format=="str"  or format=="text":
3330             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3331          else:
3332             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3333    
3334       def substitute(self,argvals):
3335          """
3336          assigns new values to symbols in the definition of the symbol.
3337          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3338    
3339          @param argvals: new values assigned to symbols
3340          @type argvals: C{dict} with keywords of type L{Symbol}.
3341          @return: result of the substitution process. Operations are executed as much as possible.
3342          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3343          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3344          """
3345          if argvals.has_key(self):
3346             arg=argvals[self]
3347             if self.isAppropriateValue(arg):
3348                return arg
3349             else:
3350                raise TypeError,"%s: new value is not appropriate."%str(self)
3351          else:
3352             arg=self.getSubstitutedArguments(argvals)
3353             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3354    
3355       def diff(self,arg):
3356          """
3357          differential of this object
3358    
3359          @param arg: the derivative is calculated with respect to arg
3360          @type arg: L{escript.Symbol}
3361          @return: derivative with respect to C{arg}
3362          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3363          """
3364          if arg==self:
3365             return identity(self.getShape())
3366          else:
3367             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3368    
3369  def symmetric(arg):  def symmetric(arg):
3370      """      """
3371      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 3173  def symmetric(arg): Line 3378  def symmetric(arg):
3378      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3379        if arg.rank==2:        if arg.rank==2:
3380          if not (arg.shape[0]==arg.shape[1]):          if not (arg.shape[0]==arg.shape[1]):
3381             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3382        elif arg.rank==4:        elif arg.rank==4:
3383          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3384             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3385        else:        else:
3386          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3387        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3388      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3389        return escript_symmetric(arg)        if arg.getRank()==2:
3390            if not (arg.getShape()[0]==arg.getShape()[1]):
3391               raise ValueError,"argument must be square."
3392            return arg._symmetric()
3393          elif arg.getRank()==4:
3394            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3395               raise ValueError,"argument must be square."
3396            return arg._symmetric()
3397          else:
3398            raise ValueError,"rank 2 or 4 is required."
3399      elif isinstance(arg,float):      elif isinstance(arg,float):
3400        return arg        return arg
3401      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3189  def symmetric(arg): Line 3403  def symmetric(arg):
3403      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3404        if arg.getRank()==2:        if arg.getRank()==2:
3405          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3406             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3407        elif arg.getRank()==4:        elif arg.getRank()==4:
3408          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3409             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3410        else:        else:
3411          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3412        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3413      else:      else:
3414        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3415    
 def escript_symmetric(arg):  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=arg._symmetric()  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=arg._symmetric()  
       else:  
         raise ValueError,"escript_symmetric: rank 2 or 4 is required."  
       return out  
   
3416  def nonsymmetric(arg):  def nonsymmetric(arg):
3417      """      """
3418      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
# Line 3232  def nonsymmetric(arg): Line 3433  def nonsymmetric(arg):
3433          raise ValueError,"nonsymmetric: rank 2 or 4 is required."          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3434        return (arg-transpose(arg))/2        return (arg-transpose(arg))/2
3435      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3436        return escript_nonsymmetric(arg)        if arg.getRank()==2:
3437            if not (arg.getShape()[0]==arg.getShape()[1]):
3438               raise ValueError,"argument must be square."
3439            return arg._nonsymmetric()
3440          elif arg.getRank()==4:
3441            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3442               raise ValueError,"argument must be square."
3443            return arg._nonsymmetric()
3444          else:
3445            raise ValueError,"rank 2 or 4 is required."
3446      elif isinstance(arg,float):      elif isinstance(arg,float):
3447        return arg        return arg
3448      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3250  def nonsymmetric(arg): Line 3460  def nonsymmetric(arg):
3460      else:      else:
3461        raise TypeError,"nonsymmetric: Unknown argument type."        raise TypeError,"nonsymmetric: Unknown argument type."
3462    
 def escript_nonsymmetric(arg): # this should be implemented in c++  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=arg._nonsymmetric()  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=arg._nonsymmetric()  
       else:  
         raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."  
       return out  
   
   
3463  def inverse(arg):  def inverse(arg):
3464      """      """
3465      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3466    
3467      @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.
3468      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3469      @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])
3470      @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
3471      @note: for L{escript.Data} objects the dimension is restricted to 3.      @note: for L{escript.Data} objects the dimension is restricted to 3.
3472      """      """
# Line 3407  class Inverse_Symbol(DependendSymbol): Line 3603  class Inverse_Symbol(DependendSymbol):
3603        if arg==self:        if arg==self:
3604           return identity(self.getShape())           return identity(self.getShape())
3605        else:        else:
3606           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3607    
3608  def eigenvalues(arg):  def eigenvalues(arg):
3609      """      """
# Line 3545  class Add_Symbol(DependendSymbol): Line 3741  class Add_Symbol(DependendSymbol):
3741         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3742         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3743         """         """
3744         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3745         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3746         if not sh0==sh1:         if not sh0==sh1:
3747            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3748         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 3620  def mult(arg0,arg1): Line 3816  def mult(arg0,arg1):
3816         """         """
3817         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3818         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3819            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3820         else:         else:
3821            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3822                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3644  class Mult_Symbol(DependendSymbol): Line 3840  class Mult_Symbol(DependendSymbol):
3840         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3841         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3842         """         """
3843         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3844         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3845         if not sh0==sh1:         if not sh0==sh1:
3846            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3847         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 3720  def quotient(arg0,arg1): Line 3916  def quotient(arg0,arg1):
3916         """         """
3917         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3918         if testForZero(args[0]):         if testForZero(args[0]):
3919            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3920         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3921            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3922               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3749  class Quotient_Symbol(DependendSymbol): Line 3945  class Quotient_Symbol(DependendSymbol):
3945         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3946         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3947         """         """
3948         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3949         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3950         if not sh0==sh1:         if not sh0==sh1:
3951            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3952         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 3826  def power(arg0,arg1): Line 4022  def power(arg0,arg1):
4022         """         """
4023         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4024         if testForZero(args[0]):         if testForZero(args[0]):
4025            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4026         elif testForZero(args[1]):         elif testForZero(args[1]):
4027            return numarray.ones(pokeShape(args[1]),numarray.Float64)            return numarray.ones(getShape(args[1]),numarray.Float64)
4028         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4029            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4030         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 3851  class Power_Symbol(DependendSymbol): Line 4047  class Power_Symbol(DependendSymbol):
4047         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4048         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4049         """         """
4050         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4051         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4052         if not sh0==sh1:         if not sh0==sh1:
4053            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4054         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3951  def minimum(*args): Line 4147  def minimum(*args):
4147            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4148      return out      return out
4149    
4150  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4151      """      """
4152      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4153    
4154      @param arg: argument      @param arg: argument
4155      @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}
4156      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4157      @type minval: C{float}      @type minval: C{float} or C{None}
4158      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4159      @type maxval: C{float}      @type maxval: C{float} or C{None}
4160      @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
4161               less then maxval are unchanged.               less then maxval are unchanged.
4162      @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
4163      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4164      """      """
4165      if minval>maxval:      if not minval==None and not maxval==None:
4166         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4167      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4168        if minval == None:
4169            tmp=arg
4170        else:
4171            tmp=maximum(minval,arg)
4172        if maxval == None:
4173            return tmp
4174        else:
4175            return minimum(tmp,maxval)
4176    
4177        
4178  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3989  def inner(arg0,arg1): Line 4193  def inner(arg0,arg1):
4193      @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
4194      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4195      """      """
4196      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4197      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4198      if not sh0==sh1:      if not sh0==sh1:
4199          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4200      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4201    
4202    def outer(arg0,arg1):
4203        """
4204        the outer product of the two argument:
4205    
4206        out[t,s]=arg0[t]*arg1[s]
4207    
4208        where
4209    
4210            - s runs through arg0.Shape
4211            - t runs through arg1.Shape
4212    
4213        @param arg0: first argument
4214        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4215        @param arg1: second argument
4216        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4217        @return: the outer product of arg0 and arg1 at each data point
4218        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4219        """
4220        return generalTensorProduct(arg0,arg1,axis_offset=0)
4221    
4222  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4223      """      """
4224        see L{matrix_mult}
4225        """
4226        return matrix_mult(arg0,arg1)
4227    
4228    def matrix_mult(arg0,arg1):
4229        """
4230      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4231    
4232      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
# Line 4005  def matrixmult(arg0,arg1): Line 4235  def matrixmult(arg0,arg1):
4235    
4236      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4237    
4238      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.
4239    
4240      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4241      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4015  def matrixmult(arg0,arg1): Line 4245  def matrixmult(arg0,arg1):
4245      @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
4246      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4247      """      """
4248      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4249      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4250      if not len(sh0)==2 :      if not len(sh0)==2 :
4251          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4252      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4253          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4254      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4255    
4256  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4257      """      """
4258      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  
4259      """      """
4260      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
4261    
4262    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4263      """      """
4264      the tensor product of the two argument:      the tensor product of the two argument:
4265            
# Line 4069  def tensormult(arg0,arg1): Line 4284  def tensormult(arg0,arg1):
4284    
4285      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]
4286    
4287      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  
4288      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.
4289    
4290      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4291      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4079  def tensormult(arg0,arg1): Line 4294  def tensormult(arg0,arg1):
4294      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4295      @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
4296      """      """
4297      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4298      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4299      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4300         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4301      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):
4302         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4303      else:      else:
4304          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4305    
4306  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4307      """      """
# Line 4100  def generalTensorProduct(arg0,arg1,axis_ Line 4315  def generalTensorProduct(arg0,arg1,axis_
4315          - r runs trough arg0.Shape[:axis_offset]          - r runs trough arg0.Shape[:axis_offset]
4316          - t runs through arg1.Shape[axis_offset:]          - t runs through arg1.Shape[axis_offset:]
4317    
     In the first case the the second dimension of arg0 and the length of arg1 must match and    
     in the second case the two last dimensions of arg0 must match the shape of arg1.  
   
4318      @param arg0: first argument      @param arg0: first argument
4319      @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}
4320      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4321      @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}
4322      @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.
4323      @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 4118  def generalTensorProduct(arg0,arg1,axis_ Line 4330  def generalTensorProduct(arg0,arg1,axis_
4330             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4331         else:         else:
4332             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4333                 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)
4334             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4335             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4336             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 4144  def generalTensorProduct(arg0,arg1,axis_ Line 4356  def generalTensorProduct(arg0,arg1,axis_
4356                                    
4357  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4358     """     """
4359     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4360     """     """
4361     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4362         """         """
4363         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4364    
4365         @param arg0: numerator         @param arg0: first argument
4366         @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}.
4367         @param arg1: denominator         @param arg1: second argument
4368         @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}.
4369         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4370         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4371         """         """
4372         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4373         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4374         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4375         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4376         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 4205  class GeneralTensorProduct_Symbol(Depend Line 4417  class GeneralTensorProduct_Symbol(Depend
4417           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4418           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4419    
4420  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4421      "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!!!"
4422      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4423      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
4424      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  def transposed_matrix_mult(arg0,arg1):
4425      shape10=arg1.getShape()[:axis_offset]      """
4426      shape1=arg1.getShape()[axis_offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4427      if not shape01==shape10:  
4428          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]
4429    
4430      # whatr function space should be used? (this here is not good!)      or
4431      fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
4432      # create return value:      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4433      out=escript.Data(0.,tuple(shape0+shape1),fs)  
4434      #      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4435      s0=[[]]  
4436      for k in shape0:      The first dimension of arg0 and arg1 must match.
4437            s=[]  
4438            for j in s0:      @param arg0: first argument of rank 2
4439                  for i in range(k): s.append(j+[slice(i,i)])      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4440            s0=s      @param arg1: second argument of at least rank 1
4441      s1=[[]]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4442      for k in shape1:      @return: the product of the transposed of arg0 and arg1 at each data point
4443            s=[]      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4444            for j in s1:      @raise ValueError: if the shapes of the arguments are not appropriate
4445                  for i in range(k): s.append(j+[slice(i,i)])      """
4446            s1=s      sh0=getShape(arg0)
4447      s01=[[]]      sh1=getShape(arg1)
4448      for k in shape01:      if not len(sh0)==2 :
4449            s=[]          raise ValueError,"first argument must have rank 2"
4450            for j in s01:      if not len(sh1)==2 and not len(sh1)==1:
4451                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"second argument must have rank 1 or 2"
4452            s01=s      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4453    
4454      for i0 in s0:  def transposed_tensor_mult(arg0,arg1):
4455         for i1 in s1:      """
4456           s=escript.Scalar(0.,fs)      the tensor product of the transposed of the first and the second argument
4457           for i01 in s01:      
4458              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      for arg0 of rank 2 this is
4459           out.__setitem__(tuple(i0+i1),s)  
4460      return out      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4461    
4462        or
4463    
4464        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4465    
4466      
4467        and for arg0 of rank 4 this is
4468    
4469        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4470                  
4471        or
4472    
4473        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4474    
4475        or
4476    
4477        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4478    
4479        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4480        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4481    
4482        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4483    
4484        @param arg0: first argument of rank 2 or 4
4485        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4486        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4487        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4488        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4489        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4490        """
4491        sh0=getShape(arg0)
4492        sh1=getShape(arg1)
4493        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4494           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4495        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4496           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4497        else:
4498            raise ValueError,"first argument must have rank 2 or 4"
4499    
4500    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4501        """
4502        generalized tensor product of transposed of arg0 and arg1:
4503    
4504        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4505    
4506        where
4507    
4508            - s runs through arg0.Shape[axis_offset:]
4509            - r runs trough arg0.Shape[:axis_offset]
4510            - t runs through arg1.Shape[axis_offset:]
4511    
4512        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4513        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4514    
4515        @param arg0: first argument
4516        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4517        @param arg1: second argument
4518        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4519        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4520        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4521        """
4522        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4523        arg0,arg1=matchType(arg0,arg1)
4524        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4525        if isinstance(arg0,numarray.NumArray):
4526           if isinstance(arg1,Symbol):
4527               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4528           else:
4529               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4530                   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)
4531               arg0_c=arg0.copy()
4532               arg1_c=arg1.copy()
4533               sh0,sh1=arg0.shape,arg1.shape
4534               d0,d1,d01=1,1,1
4535               for i in sh0[axis_offset:]: d0*=i
4536               for i in sh1[axis_offset:]: d1*=i
4537               for i in sh0[:axis_offset]: d01*=i
4538               arg0_c.resize((d01,d0))
4539               arg1_c.resize((d01,d1))
4540               out=numarray.zeros((d0,d1),numarray.Float64)
4541               for i0 in range(d0):
4542                        for i1 in range(d1):
4543                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4544               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4545               return out
4546        elif isinstance(arg0,escript.Data):
4547           if isinstance(arg1,Symbol):
4548               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4549           else:
4550               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4551        else:      
4552           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4553                    
4554    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4555       """
4556       Symbol representing the general tensor product of the transposed of arg0 and arg1
4557       """
4558       def __init__(self,arg0,arg1,axis_offset=0):
4559           """
4560           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4561    
4562           @param arg0: first argument
4563           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4564           @param arg1: second argument
4565           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4566           @raise ValueError: inconsistent dimensions of arguments.
4567           @note: if both arguments have a spatial dimension, they must equal.
4568           """
4569           sh_arg0=getShape(arg0)
4570           sh_arg1=getShape(arg1)
4571           sh01=sh_arg0[:axis_offset]
4572           sh10=sh_arg1[:axis_offset]
4573           sh0=sh_arg0[axis_offset:]
4574           sh1=sh_arg1[axis_offset:]
4575           if not sh01==sh10:
4576               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)
4577           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4578    
4579       def getMyCode(self,argstrs,format="escript"):
4580          """
4581          returns a program code that can be used to evaluate the symbol.
4582    
4583          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4584          @type argstrs: C{list} of length 2 of C{str}.
4585          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4586          @type format: C{str}
4587          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4588          @rtype: C{str}
4589          @raise NotImplementedError: if the requested format is not available
4590          """
4591          if format=="escript" or format=="str" or format=="text":
4592             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4593          else:
4594             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4595    
4596       def substitute(self,argvals):
4597          """
4598          assigns new values to symbols in the definition of the symbol.
4599          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4600    
4601          @param argvals: new values assigned to symbols
4602          @type argvals: C{dict} with keywords of type L{Symbol}.
4603          @return: result of the substitution process. Operations are executed as much as possible.
4604          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4605          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4606          """
4607          if argvals.has_key(self):
4608             arg=argvals[self]
4609             if self.isAppropriateValue(arg):
4610                return arg
4611             else:
4612                raise TypeError,"%s: new value is not appropriate."%str(self)
4613          else:
4614             args=self.getSubstitutedArguments(argvals)
4615             return generalTransposedTensorProduct(args[0],args[1],args[2])
4616    
4617    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4618        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4619        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4620    
4621    def matrix_transposed_mult(arg0,arg1):
4622        """
4623        matrix-transposed(matrix) product of the two argument:
4624    
4625        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4626    
4627        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4628    
4629        The last dimensions of arg0 and arg1 must match.
4630    
4631        @param arg0: first argument of rank 2
4632        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4633        @param arg1: second argument of rank 2
4634        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4635        @return: the product of arg0 and the transposed of arg1 at each data point
4636        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4637        @raise ValueError: if the shapes of the arguments are not appropriate
4638        """
4639        sh0=getShape(arg0)
4640        sh1=getShape(arg1)
4641        if not len(sh0)==2 :
4642            raise ValueError,"first argument must have rank 2"
4643        if not len(sh1)==2 and not len(sh1)==1:
4644            raise ValueError,"second argument must have rank 1 or 2"
4645        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4646    
4647    def tensor_transposed_mult(arg0,arg1):
4648        """
4649        the tensor product of the first and the transpose of the second argument
4650        
4651        for arg0 of rank 2 this is
4652    
4653        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4654    
4655        and for arg0 of rank 4 this is
4656    
4657        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4658                  
4659        or
4660    
4661        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4662    
4663        In the first case the the second dimension of arg0 and arg1 must match and  
4664        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4665    
4666        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4667    
4668        @param arg0: first argument of rank 2 or 4
4669        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4670        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4671        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4672        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4673        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4674        """
4675        sh0=getShape(arg0)
4676        sh1=getShape(arg1)
4677        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4678           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4679        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4680           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4681        else:
4682            raise ValueError,"first argument must have rank 2 or 4"
4683    
4684    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4685        """
4686        generalized tensor product of transposed of arg0 and arg1:
4687    
4688        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4689    
4690        where
4691    
4692            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4693            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4694            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4695    
4696        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4697        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4698    
4699        @param arg0: first argument
4700        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4701        @param arg1: second argument
4702        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4703        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4704        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4705        """
4706        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4707        arg0,arg1=matchType(arg0,arg1)
4708        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4709        if isinstance(arg0,numarray.NumArray):
4710           if isinstance(arg1,Symbol):
4711               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4712           else:
4713               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4714                   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)
4715               arg0_c=arg0.copy()
4716               arg1_c=arg1.copy()
4717               sh0,sh1=arg0.shape,arg1.shape
4718               d0,d1,d01=1,1,1
4719               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4720               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4721               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4722               arg0_c.resize((d0,d01))
4723               arg1_c.resize((d1,d01))
4724               out=numarray.zeros((d0,d1),numarray.Float64)
4725               for i0 in range(d0):
4726                        for i1 in range(d1):
4727                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4728               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4729               return out
4730        elif isinstance(arg0,escript.Data):
4731           if isinstance(arg1,Symbol):
4732               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4733           else:
4734               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4735        else:      
4736           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4737                    
4738    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4739       """
4740       Symbol representing the general tensor product of arg0 and the transpose of arg1
4741       """
4742       def __init__(self,arg0,arg1,axis_offset=0):
4743           """
4744           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4745    
4746           @param arg0: first argument
4747           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4748           @param arg1: second argument
4749           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4750           @raise ValueError: inconsistent dimensions of arguments.
4751           @note: if both arguments have a spatial dimension, they must equal.
4752           """
4753           sh_arg0=getShape(arg0)
4754           sh_arg1=getShape(arg1)
4755           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4756           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4757           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4758           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4759           if not sh01==sh10:
4760               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)
4761           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4762    
4763       def getMyCode(self,argstrs,format="escript"):
4764          """
4765          returns a program code that can be used to evaluate the symbol.
4766    
4767          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4768          @type argstrs: C{list} of length 2 of C{str}.
4769          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4770          @type format: C{str}
4771          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4772          @rtype: C{str}
4773          @raise NotImplementedError: if the requested format is not available
4774          """
4775          if format=="escript" or format=="str" or format=="text":
4776             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4777          else:
4778             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4779    
4780       def substitute(self,argvals):
4781          """
4782          assigns new values to symbols in the definition of the symbol.
4783          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4784    
4785          @param argvals: new values assigned to symbols
4786          @type argvals: C{dict} with keywords of type L{Symbol}.
4787          @return: result of the substitution process. Operations are executed as much as possible.
4788          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4789          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4790          """
4791          if argvals.has_key(self):
4792             arg=argvals[self]
4793             if self.isAppropriateValue(arg):
4794                return arg
4795             else:
4796                raise TypeError,"%s: new value is not appropriate."%str(self)
4797          else:
4798             args=self.getSubstitutedArguments(argvals)
4799             return generalTensorTransposedProduct(args[0],args[1],args[2])
4800    
4801    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4802        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4803        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4804    
4805  #=========================================================  #=========================================================
4806  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
# Line 4372  def integrate(arg,where=None): Line 4926  def integrate(arg,where=None):
4926         else:         else:
4927            return arg._integrate()            return arg._integrate()
4928      else:      else:
4929        raise TypeError,"integrate: Unknown argument type."         arg2=escript.Data(arg,where)
4930           if arg2.getRank()==0:
4931              return arg2._integrate()[0]
4932           else:
4933              return arg2._integrate()
4934    
4935  class Integrate_Symbol(DependendSymbol):  class Integrate_Symbol(DependendSymbol):
4936     """     """

Legend:
Removed from v.775  
changed lines
  Added in v.1295

  ViewVC Help
Powered by ViewVC 1.1.26