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

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

Parent Directory | Revision Log | View Patch Patch

-revision 775 by ksteube,
Mon Jul 10 04:00:08 2006 UTC
+revision 1044 by gross,
Mon Mar 19 07:29:31 2007 UTC
 Line 26 
 import math
  import numarray
  import escript
  import os
+ from esys.escript import C_GeneralTensorProduct
  #=========================================================
  #   some helpers:
  #=========================================================
+ def getTagNames(domain):
+     """
+     returns a list of the tag names used by the domain
+     @param domain: a domain object
+     @type domain: L{escript.Domain}
+     @return: a list of the tag name used by the domain.
+     @rtype: C{list} of C{str}
+     """
+     return [n.strip() for n in domain.showTagNames().split(",") ]
+ def insertTagNames(domain,**kwargs):
+     """
+     inserts tag names into the domain
+     @param domain: a domain object
+     @type domain: C{escript.Domain}
+     @keyword <tag name>: tag key assigned to <tag name>
+     @type <tag name>: C{int}
+     """
+     for  k in kwargs:
+          domain.setTagMap(k,kwargs[k])
+ def insertTaggedValues(target,**kwargs):
+     """
+     inserts tagged values into the tagged using tag names
+     @param target: data to be filled by tagged values
+     @type target: L{escript.Data}
+     @keyword <tag name>: value to be used for <tag name>
+     @type <tag name>: C{float} or {numarray.NumArray}
+     @return: C{target}
+     @rtype: L{escript.Data}
+     """
+     for k in kwargs:
+         target.setTaggedValue(k,kwargs[k])
+     return target
  def saveVTK(filename,domain=None,**data):
      """
      writes a L{Data} objects into a files using the the VTK XML file format.
-Line 238 
 def inf(arg):
+Line 279 
 def inf(arg):
  #=========================================================================
  #   some little helpers
  #=========================================================================
- def pokeShape(arg):
+ def getRank(arg):
+     """
+     identifies the rank of its argument
+     @param arg: a given object
+     @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
+     @return: the rank of the argument
+     @rtype: C{int}
+     @raise TypeError: if type of arg cannot be processed
+     """
+     if isinstance(arg,numarray.NumArray):
+         return arg.rank
+     elif isinstance(arg,escript.Data):
+         return arg.getRank()
+     elif isinstance(arg,float):
+         return 0
+     elif isinstance(arg,int):
+         return 0
+     elif isinstance(arg,Symbol):
+         return arg.getRank()
+     else:
+       raise TypeError,"getShape: cannot identify shape"
+ def getShape(arg):
      """
      identifies the shape of its argument
-Line 260 
 def pokeShape(arg):
+Line 324 
 def pokeShape(arg):
      elif isinstance(arg,Symbol):
          return arg.getShape()
      else:
-       raise TypeError,"pokeShape: cannot identify shape"
+       raise TypeError,"getShape: cannot identify shape"
  def pokeDim(arg):
      """
-Line 283 
 def commonShape(arg0,arg1):
+Line 347 
 def commonShape(arg0,arg1):
      """
      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.
-     @param arg0: an object with a shape (see L{pokeShape})
+     @param arg0: an object with a shape (see L{getShape})
-     @param arg1: an object with a shape (see L{pokeShape})
+     @param arg1: an object with a shape (see L{getShape})
      @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.
      @rtype: C{tuple} of C{int}
      @raise ValueError: if no shape can be found.
      """
-     sh0=pokeShape(arg0)
+     sh0=getShape(arg0)
-     sh1=pokeShape(arg1)
+     sh1=getShape(arg1)
      if len(sh0)<len(sh1):
         if not sh0==sh1[:len(sh0)]:
               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
-Line 440 
 def matchShape(arg0,arg1):
+Line 504 
 def matchShape(arg0,arg1):
      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
      @param arg1: a given object
      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
-     @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.
      @rtype: C{tuple}
      """
      sh=commonShape(arg0,arg1)
-     sh0=pokeShape(arg0)
+     sh0=getShape(arg0)
-     sh1=pokeShape(arg1)
+     sh1=getShape(arg1)
      if len(sh0)<len(sh):
         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
      elif len(sh1)<len(sh):
-Line 573 
 class Symbol(object):
+Line 637 
 class Symbol(object):
            if isinstance(a,Symbol):
               out.append(a.substitute(argvals))
            else:
-               s=pokeShape(s)+arg.getShape()
+               s=getShape(s)+arg.getShape()
                if len(s)>0:
                   out.append(numarray.zeros(s),numarray.Float64)
                else:
-Line 1310 
 def whereNonZero(arg,tol=0.):
+Line 1374 
 def whereNonZero(arg,tol=0.):
     else:
        raise TypeError,"whereNonZero: Unknown argument type."
+ def erf(arg):
+    """
+    returns erf of argument arg
+    @param arg: argument
+    @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
+    @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
+    @raises TypeError: if the type of the argument is not expected.
+    """
+    if isinstance(arg,escript.Data):
+       return arg._erf()
+    else:
+       raise TypeError,"erf: Unknown argument type."
  def sin(arg):
     """
     returns sine of argument arg
-Line 2930 
 def trace(arg,axis_offset=0):
+Line 3008 
 def trace(arg,axis_offset=0):
     @param arg: argument
     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
-    @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
-                   axis_offset and axis_offset+1 must be equal.
+                   C{axis_offset} and axis_offset+1 must be equal.
     @type axis_offset: C{int}
     @return: trace of arg. The rank of the returned object is minus 2 of the rank 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 3017 
 def trace(arg,axis_offset=0):
     if isinstance(arg,numarray.NumArray):
        sh=arg.shape
        if len(sh)<2:
-         raise ValueError,"trace: rank of argument must be greater than 1"
+         raise ValueError,"rank of argument must be greater than 1"
        if axis_offset<0 or axis_offset>len(sh)-2:
-         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
        s1=1
        for i in range(axis_offset): s1*=sh[i]
        s2=1
        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
        if not sh[axis_offset] == sh[axis_offset+1]:
-         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)
        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
        out=numarray.zeros([s1,s2],numarray.Float64)
        for i1 in range(s1):
-Line 2956 
 def trace(arg,axis_offset=0):
+Line 3034 
 def trace(arg,axis_offset=0):
        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
        return out
     elif isinstance(arg,escript.Data):
-       return escript_trace(arg,axis_offset)
+       if arg.getRank()<2:
+         raise ValueError,"rank of argument must be greater than 1"
+       if axis_offset<0 or axis_offset>arg.getRank()-2:
+         raise ValueError,"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,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
+       return arg._trace(axis_offset)
     elif isinstance(arg,float):
-       raise TypeError,"trace: illegal argument type float."
+       raise TypeError,"illegal argument type float."
     elif isinstance(arg,int):
-       raise TypeError,"trace: illegal argument type int."
+       raise TypeError,"illegal argument type int."
     elif isinstance(arg,Symbol):
        return Trace_Symbol(arg,axis_offset)
     else:
-       raise TypeError,"trace: Unknown argument type."
+       raise TypeError,"Unknown argument type."
- 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
  class Trace_Symbol(DependendSymbol):
     """
     L{Symbol} representing the result of the trace function
-Line 2986 
 class Trace_Symbol(DependendSymbol):
+Line 3060 
 class Trace_Symbol(DependendSymbol):
        initialization of trace L{Symbol} with argument arg
        @param arg: argument of function
        @type arg: L{Symbol}.
-       @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
-                   axis_offset and axis_offset+1 must be equal.
+                   C{axis_offset} and axis_offset+1 must be equal.
        @type axis_offset: C{int}
        """
        if arg.getRank()<2:
-         raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"
+         raise ValueError,"rank of argument must be greater than 1"
        if axis_offset<0 or axis_offset>arg.getRank()-2:
-         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
        s=list(arg.getShape())
        if not s[axis_offset] == s[axis_offset+1]:
-         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)
        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
     def getMyCode(self,argstrs,format="escript"):
-Line 3053 
 class Trace_Symbol(DependendSymbol):
+Line 3127 
 class Trace_Symbol(DependendSymbol):
  def transpose(arg,axis_offset=None):
     """
-    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.
     @param arg: argument
     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
-    @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.
-                        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.
     @type axis_offset: C{int}
     @return: transpose 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 3141 
 def transpose(arg,axis_offset=None):
        if axis_offset==None: axis_offset=int(arg.rank/2)
        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
     elif isinstance(arg,escript.Data):
-       if axis_offset==None: axis_offset=int(arg.getRank()/2)
+       r=arg.getRank()
-       return escript_transpose(arg,axis_offset)
+       if axis_offset==None: axis_offset=int(r/2)
+       if axis_offset<0 or axis_offset>r:
+         raise ValueError,"axis_offset must be between 0 and %s"%r
+       return arg._transpose(axis_offset)
     elif isinstance(arg,float):
        if not ( axis_offset==0 or axis_offset==None):
-         raise ValueError,"transpose: axis_offset must be 0 for float argument"
+         raise ValueError,"axis_offset must be 0 for float argument"
        return arg
     elif isinstance(arg,int):
        if not ( axis_offset==0 or axis_offset==None):
-         raise ValueError,"transpose: axis_offset must be 0 for int argument"
+         raise ValueError,"axis_offset must be 0 for int argument"
        return float(arg)
     elif isinstance(arg,Symbol):
        if axis_offset==None: axis_offset=int(arg.getRank()/2)
        return Transpose_Symbol(arg,axis_offset)
     else:
-       raise TypeError,"transpose: Unknown argument type."
+       raise TypeError,"Unknown argument type."
- 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
  class Transpose_Symbol(DependendSymbol):
     """
     L{Symbol} representing the result of the transpose function
-Line 3100 
 class Transpose_Symbol(DependendSymbol):
+Line 3170 
 class Transpose_Symbol(DependendSymbol):
        @param arg: argument of function
        @type arg: L{Symbol}.
-       @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.
-                        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.
        @type axis_offset: C{int}
        """
        if axis_offset==None: axis_offset=int(arg.getRank()/2)
        if axis_offset<0 or axis_offset>arg.getRank():
-         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()
        s=arg.getShape()
        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
-Line 3161 
 class Transpose_Symbol(DependendSymbol):
+Line 3231 
 class Transpose_Symbol(DependendSymbol):
           return identity(self.getShape())
        else:
           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
+ def swap_axes(arg,axis0=0,axis1=1):
+    """
+    returns the swap of arg by swaping the components axis0 and axis1
+    @param arg: argument
+    @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
+    @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
+    @type axis0: C{int}
+    @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
+    @type axis1: C{int}
+    @return: C{arg} with swaped components
+    @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
+    """
+    if axis0 > axis1:
+       axis0,axis1=axis1,axis0
+    if isinstance(arg,numarray.NumArray):
+       return numarray.swapaxes(arg,axis0,axis1)
+    elif isinstance(arg,escript.Data):
+       return arg._swap_axes(axis0,axis1)
+    elif isinstance(arg,float):
+       raise TyepError,"float argument is not supported."
+    elif isinstance(arg,int):
+       raise TyepError,"int argument is not supported."
+    elif isinstance(arg,Symbol):
+       return SwapAxes_Symbol(arg,axis0,axis1)
+    else:
+       raise TypeError,"Unknown argument type."
+ class SwapAxes_Symbol(DependendSymbol):
+    """
+    L{Symbol} representing the result of the swap function
+    """
+    def __init__(self,arg,axis0=0,axis1=1):
+       """
+       initialization of swap L{Symbol} with argument arg
+       @param arg: argument
+       @type arg: L{Symbol}.
+       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
+       @type axis0: C{int}
+       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
+       @type axis1: C{int}
+       """
+       if arg.getRank()<2:
+          raise ValueError,"argument must have at least rank 2."
+       if axis0<0 or axis0>arg.getRank()-1:
+          raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
+       if axis1<0 or axis1>arg.getRank()-1:
+          raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
+       if axis0 == axis1:
+          raise ValueError,"axis indices must be different."
+       if axis0 > axis1:
+          axis0,axis1=axis1,axis0
+       s=arg.getShape()
+       s_out=[]
+       for i in range(len(s)):
+          if i == axis0:
+             s_out.append(s[axis1])
+          elif i == axis1:
+             s_out.append(s[axis0])
+          else:
+             s_out.append(s[i])
+       super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
+    def getMyCode(self,argstrs,format="escript"):
+       """
+       returns a program code that can be used to evaluate the symbol.
+       @param argstrs: gives for each argument a string representing the argument for the evaluation.
+       @type argstrs: C{str} or a C{list} of length 1 of C{str}.
+       @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
+       @type format: C{str}
+       @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
+       @rtype: C{str}
+       @raise NotImplementedError: if the requested format is not available
+       """
+       if format=="escript" or format=="str"  or format=="text":
+          return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
+       else:
+          raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
+    def substitute(self,argvals):
+       """
+       assigns new values to symbols in the definition of the symbol.
+       The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
+       @param argvals: new values assigned to symbols
+       @type argvals: C{dict} with keywords of type L{Symbol}.
+       @return: result of the substitution process. Operations are executed as much as possible.
+       @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
+       @raise TypeError: if a value for a L{Symbol} cannot be substituted.
+       """
+       if argvals.has_key(self):
+          arg=argvals[self]
+          if self.isAppropriateValue(arg):
+             return arg
+          else:
+             raise TypeError,"%s: new value is not appropriate."%str(self)
+       else:
+          arg=self.getSubstitutedArguments(argvals)
+          return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
+    def diff(self,arg):
+       """
+       differential of this object
+       @param arg: the derivative is calculated with respect to arg
+       @type arg: L{escript.Symbol}
+       @return: derivative with respect to C{arg}
+       @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
+       """
+       if arg==self:
+          return identity(self.getShape())
+       else:
+          return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
  def symmetric(arg):
      """
      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
-Line 3173 
 def symmetric(arg):
+Line 3360 
 def symmetric(arg):
      if isinstance(arg,numarray.NumArray):
        if arg.rank==2:
          if not (arg.shape[0]==arg.shape[1]):
-            raise ValueError,"symmetric: argument must be square."
+            raise ValueError,"argument must be square."
        elif arg.rank==4:
          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
-            raise ValueError,"symmetric: argument must be square."
+            raise ValueError,"argument must be square."
        else:
-         raise ValueError,"symmetric: rank 2 or 4 is required."
+         raise ValueError,"rank 2 or 4 is required."
        return (arg+transpose(arg))/2
      elif isinstance(arg,escript.Data):
-       return escript_symmetric(arg)
+       if arg.getRank()==2:
+         if not (arg.getShape()[0]==arg.getShape()[1]):
+            raise ValueError,"argument must be square."
+         return arg._symmetric()
+       elif arg.getRank()==4:
+         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
+            raise ValueError,"argument must be square."
+         return arg._symmetric()
+       else:
+         raise ValueError,"rank 2 or 4 is required."
      elif isinstance(arg,float):
        return arg
      elif isinstance(arg,int):
-Line 3189 
 def symmetric(arg):
+Line 3385 
 def symmetric(arg):
      elif isinstance(arg,Symbol):
        if arg.getRank()==2:
          if not (arg.getShape()[0]==arg.getShape()[1]):
-            raise ValueError,"symmetric: argument must be square."
+            raise ValueError,"argument must be square."
        elif arg.getRank()==4:
          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
-            raise ValueError,"symmetric: argument must be square."
+            raise ValueError,"argument must be square."
        else:
-         raise ValueError,"symmetric: rank 2 or 4 is required."
+         raise ValueError,"rank 2 or 4 is required."
        return (arg+transpose(arg))/2
      else:
        raise TypeError,"symmetric: Unknown argument type."
- 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
  def nonsymmetric(arg):
      """
      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
-Line 3232 
 def nonsymmetric(arg):
+Line 3415 
 def nonsymmetric(arg):
          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
        return (arg-transpose(arg))/2
      elif isinstance(arg,escript.Data):
-       return escript_nonsymmetric(arg)
+       if arg.getRank()==2:
+         if not (arg.getShape()[0]==arg.getShape()[1]):
+            raise ValueError,"argument must be square."
+         return arg._nonsymmetric()
+       elif arg.getRank()==4:
+         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
+            raise ValueError,"argument must be square."
+         return arg._nonsymmetric()
+       else:
+         raise ValueError,"rank 2 or 4 is required."
      elif isinstance(arg,float):
        return arg
      elif isinstance(arg,int):
-Line 3250 
 def nonsymmetric(arg):
+Line 3442 
 def nonsymmetric(arg):
      else:
        raise TypeError,"nonsymmetric: Unknown argument type."
- 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
  def inverse(arg):
      """
      returns the inverse of the square matrix arg.
      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
-     @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])
      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
      @note: for L{escript.Data} objects the dimension is restricted to 3.
      """
-Line 3407 
 class Inverse_Symbol(DependendSymbol):
+Line 3585 
 class Inverse_Symbol(DependendSymbol):
        if arg==self:
           return identity(self.getShape())
        else:
-          return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)
+          return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
  def eigenvalues(arg):
      """
-Line 3545 
 class Add_Symbol(DependendSymbol):
+Line 3723 
 class Add_Symbol(DependendSymbol):
         @raise ValueError: if both arguments do not have the same shape.
         @note: if both arguments have a spatial dimension, they must equal.
         """
-        sh0=pokeShape(arg0)
+        sh0=getShape(arg0)
-        sh1=pokeShape(arg1)
+        sh1=getShape(arg1)
         if not sh0==sh1:
            raise ValueError,"Add_Symbol: shape of arguments must match"
         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
-Line 3620 
 def mult(arg0,arg1):
+Line 3798 
 def mult(arg0,arg1):
         """
         args=matchShape(arg0,arg1)
         if testForZero(args[0]) or testForZero(args[1]):
-           return numarray.zeros(pokeShape(args[0]),numarray.Float64)
+           return numarray.zeros(getShape(args[0]),numarray.Float64)
         else:
            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
                return Mult_Symbol(args[0],args[1])
-Line 3644 
 class Mult_Symbol(DependendSymbol):
+Line 3822 
 class Mult_Symbol(DependendSymbol):
         @raise ValueError: if both arguments do not have the same shape.
         @note: if both arguments have a spatial dimension, they must equal.
         """
-        sh0=pokeShape(arg0)
+        sh0=getShape(arg0)
-        sh1=pokeShape(arg1)
+        sh1=getShape(arg1)
         if not sh0==sh1:
            raise ValueError,"Mult_Symbol: shape of arguments must match"
         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
-Line 3720 
 def quotient(arg0,arg1):
+Line 3898 
 def quotient(arg0,arg1):
         """
         args=matchShape(arg0,arg1)
         if testForZero(args[0]):
-           return numarray.zeros(pokeShape(args[0]),numarray.Float64)
+           return numarray.zeros(getShape(args[0]),numarray.Float64)
         elif isinstance(args[0],Symbol):
            if isinstance(args[1],Symbol):
               return Quotient_Symbol(args[0],args[1])
-Line 3749 
 class Quotient_Symbol(DependendSymbol):
+Line 3927 
 class Quotient_Symbol(DependendSymbol):
         @raise ValueError: if both arguments do not have the same shape.
         @note: if both arguments have a spatial dimension, they must equal.
         """
-        sh0=pokeShape(arg0)
+        sh0=getShape(arg0)
-        sh1=pokeShape(arg1)
+        sh1=getShape(arg1)
         if not sh0==sh1:
            raise ValueError,"Quotient_Symbol: shape of arguments must match"
         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
-Line 3826 
 def power(arg0,arg1):
+Line 4004 
 def power(arg0,arg1):
         """
         args=matchShape(arg0,arg1)
         if testForZero(args[0]):
-           return numarray.zeros(pokeShape(args[0]),numarray.Float64)
+           return numarray.zeros(getShape(args[0]),numarray.Float64)
         elif testForZero(args[1]):
-           return numarray.ones(pokeShape(args[1]),numarray.Float64)
+           return numarray.ones(getShape(args[1]),numarray.Float64)
         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
            return Power_Symbol(args[0],args[1])
         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):
-Line 3851 
 class Power_Symbol(DependendSymbol):
+Line 4029 
 class Power_Symbol(DependendSymbol):
         @raise ValueError: if both arguments do not have the same shape.
         @note: if both arguments have a spatial dimension, they must equal.
         """
-        sh0=pokeShape(arg0)
+        sh0=getShape(arg0)
-        sh1=pokeShape(arg1)
+        sh1=getShape(arg1)
         if not sh0==sh1:
            raise ValueError,"Power_Symbol: shape of arguments must match"
         d0=pokeDim(arg0)
-Line 3951 
 def minimum(*args):
+Line 4129 
 def minimum(*args):
            out=add(out,mult(whereNegative(diff),diff))
      return out
- def clip(arg,minval=0.,maxval=1.):
+ def clip(arg,minval=None,maxval=None):
      """
      cuts the values of arg between minval and maxval
      @param arg: argument
      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
-     @param minval: lower range
+     @param minval: lower range. If None no lower range is applied
-     @type minval: C{float}
+     @type minval: C{float} or C{None}
-     @param maxval: upper range
+     @param maxval: upper range. If None no upper range is applied
-     @type maxval: C{float}
+     @type maxval: C{float} or C{None}
      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
               less then maxval are unchanged.
      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
      @raise ValueError: if minval>maxval
      """
-     if minval>maxval:
+     if not minval==None and not maxval==None:
-        raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
+        if minval>maxval:
-     return minimum(maximum(minval,arg),maxval)
+           raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
+     if minval == None:
+         tmp=arg
+     else:
+         tmp=maximum(minval,arg)
+     if maxval == None:
+         return tmp
+     else:
+         return minimum(tmp,maxval)
  def inner(arg0,arg1):
-Line 3989 
 def inner(arg0,arg1):
+Line 4175 
 def inner(arg0,arg1):
      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input
      @raise ValueError: if the shapes of the arguments are not identical
      """
-     sh0=pokeShape(arg0)
+     sh0=getShape(arg0)
-     sh1=pokeShape(arg1)
+     sh1=getShape(arg1)
      if not sh0==sh1:
          raise ValueError,"inner: shape of arguments does not match"
      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
+ def outer(arg0,arg1):
+     """
+     the outer product of the two argument:
+     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
+     """
+     return generalTensorProduct(arg0,arg1,axis_offset=0)
  def matrixmult(arg0,arg1):
      """
+     see L{matrix_mult}
+     """
+     return matrix_mult(arg0,arg1)
+ def matrix_mult(arg0,arg1):
+     """
      matrix-matrix or matrix-vector product of the two argument:
      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
-Line 4005 
 def matrixmult(arg0,arg1):
+Line 4217 
 def matrixmult(arg0,arg1):
      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
-     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.
      @param arg0: first argument of rank 2
      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
-Line 4015 
 def matrixmult(arg0,arg1):
+Line 4227 
 def matrixmult(arg0,arg1):
      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
      @raise ValueError: if the shapes of the arguments are not appropriate
      """
-     sh0=pokeShape(arg0)
+     sh0=getShape(arg0)
-     sh1=pokeShape(arg1)
+     sh1=getShape(arg1)
      if not len(sh0)==2 :
          raise ValueError,"first argument must have rank 2"
      if not len(sh1)==2 and not len(sh1)==1:
          raise ValueError,"second argument must have rank 1 or 2"
      return generalTensorProduct(arg0,arg1,axis_offset=1)
- def outer(arg0,arg1):
+ def tensormult(arg0,arg1):
      """
-     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
      """
-     return generalTensorProduct(arg0,arg1,axis_offset=0)
+     return tensor_mult(arg0,arg1)
- def tensormult(arg0,arg1):
+ def tensor_mult(arg0,arg1):
      """
      the tensor product of the two argument:
-Line 4069 
 def tensormult(arg0,arg1):
+Line 4266 
 def tensormult(arg0,arg1):
      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]
-     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
-     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.
      @param arg0: first argument of rank 2 or 4
      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
-Line 4079 
 def tensormult(arg0,arg1):
+Line 4276 
 def tensormult(arg0,arg1):
      @return: the tensor product of arg0 and arg1 at each data point
      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
      """
-     sh0=pokeShape(arg0)
+     sh0=getShape(arg0)
-     sh1=pokeShape(arg1)
+     sh1=getShape(arg1)
      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
         return generalTensorProduct(arg0,arg1,axis_offset=1)
      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
         return generalTensorProduct(arg0,arg1,axis_offset=2)
      else:
-         raise ValueError,"tensormult: first argument must have rank 2 or 4"
+         raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
  def generalTensorProduct(arg0,arg1,axis_offset=0):
      """
-Line 4100 
 def generalTensorProduct(arg0,arg1,axis_
+Line 4297 
 def generalTensorProduct(arg0,arg1,axis_
          - r runs trough arg0.Shape[:axis_offset]
          - t runs through arg1.Shape[axis_offset:]
-     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.
      @param arg0: first argument
      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
-     @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
+     @param arg1: second argument
      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
      @return: the general tensor product of arg0 and arg1 at each data point.
      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
-Line 4118 
 def generalTensorProduct(arg0,arg1,axis_
+Line 4312 
 def generalTensorProduct(arg0,arg1,axis_
             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
         else:
             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
-                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)
             arg0_c=arg0.copy()
             arg1_c=arg1.copy()
             sh0,sh1=arg0.shape,arg1.shape
-Line 4144 
 def generalTensorProduct(arg0,arg1,axis_
+Line 4338 
 def generalTensorProduct(arg0,arg1,axis_
  class GeneralTensorProduct_Symbol(DependendSymbol):
     """
-    Symbol representing the quotient of two arguments.
+    Symbol representing the general tensor product of two arguments
     """
     def __init__(self,arg0,arg1,axis_offset=0):
         """
-        initialization of L{Symbol} representing the quotient of two arguments
+        initialization of L{Symbol} representing the general tensor product of two arguments.
-        @param arg0: numerator
+        @param arg0: first argument
         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @param arg1: denominator
+        @param arg1: second argument
         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
-        @raise ValueError: if both arguments do not have the same shape.
+        @raise ValueError: illegal dimension
         @note: if both arguments have a spatial dimension, they must equal.
         """
-        sh_arg0=pokeShape(arg0)
+        sh_arg0=getShape(arg0)
-        sh_arg1=pokeShape(arg1)
+        sh_arg1=getShape(arg1)
         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
         sh10=sh_arg1[:axis_offset]
-Line 4205 
 class GeneralTensorProduct_Symbol(Depend
+Line 4399 
 class GeneralTensorProduct_Symbol(Depend
           args=self.getSubstitutedArguments(argvals)
           return generalTensorProduct(args[0],args[1],args[2])
- def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct
+ def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
-     # 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:]
+ def transposed_matrix_mult(arg0,arg1):
-     shape10=arg1.getShape()[:axis_offset]
+     """
-     shape1=arg1.getShape()[axis_offset:]
+     transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
-     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)
+     out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
-     # whatr function space should be used? (this here is not good!)
+     or
-     fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
-     # create return value:
+     out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
-     out=escript.Data(0.,tuple(shape0+shape1),fs)
-     #
+     The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
-     s0=[[]]
-     for k in shape0:
+     The first dimension of arg0 and arg1 must match.
-           s=[]
-           for j in s0:
+     @param arg0: first argument of rank 2
-                 for i in range(k): s.append(j+[slice(i,i)])
+     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
-           s0=s
+     @param arg1: second argument of at least rank 1
-     s1=[[]]
+     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
-     for k in shape1:
+     @return: the product of the transposed of arg0 and arg1 at each data point
-           s=[]
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
-           for j in s1:
+     @raise ValueError: if the shapes of the arguments are not appropriate
-                 for i in range(k): s.append(j+[slice(i,i)])
+     """
-           s1=s
+     sh0=getShape(arg0)
-     s01=[[]]
+     sh1=getShape(arg1)
-     for k in shape01:
+     if not len(sh0)==2 :
-           s=[]
+         raise ValueError,"first argument must have rank 2"
-           for j in s01:
+     if not len(sh1)==2 and not len(sh1)==1:
-                 for i in range(k): s.append(j+[slice(i,i)])
+         raise ValueError,"second argument must have rank 1 or 2"
-           s01=s
+     return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
-     for i0 in s0:
+ def transposed_tensor_mult(arg0,arg1):
-        for i1 in s1:
+     """
-          s=escript.Scalar(0.,fs)
+     the tensor product of the transposed of the first and the second argument
-          for i01 in s01:
-             s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))
+     for arg0 of rank 2 this is
-          out.__setitem__(tuple(i0+i1),s)
-     return out
+     out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
+     or
+     out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
+     and for arg0 of rank 4 this is
+     out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
+     or
+     out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
+     or
+     out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
+     In the first case the the first dimension of arg0 and the first dimension of arg1 must match and
+     in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
+     The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
+     @param arg0: first argument of rank 2 or 4
+     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
+     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     """
+     sh0=getShape(arg0)
+     sh1=getShape(arg1)
+     if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
+        return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
+     elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
+        return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
+     else:
+         raise ValueError,"first argument must have rank 2 or 4"
+ def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
+     """
+     generalized tensor product of transposed of arg0 and arg1:
+     out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
+     where
+         - s runs through arg0.Shape[axis_offset:]
+         - r runs trough arg0.Shape[:axis_offset]
+         - t runs through arg1.Shape[axis_offset:]
+     The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
+     to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
+     @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 general tensor product of transposed(arg0) and arg1 at each data point.
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     """
+     if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
+     arg0,arg1=matchType(arg0,arg1)
+     # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
+     if isinstance(arg0,numarray.NumArray):
+        if isinstance(arg1,Symbol):
+            return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
+        else:
+            if not arg0.shape[:axis_offset]==arg1.shape[: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)
+            arg0_c=arg0.copy()
+            arg1_c=arg1.copy()
+            sh0,sh1=arg0.shape,arg1.shape
+            d0,d1,d01=1,1,1
+            for i in sh0[axis_offset:]: d0*=i
+            for i in sh1[axis_offset:]: d1*=i
+            for i in sh0[:axis_offset]: d01*=i
+            arg0_c.resize((d01,d0))
+            arg1_c.resize((d01,d1))
+            out=numarray.zeros((d0,d1),numarray.Float64)
+            for i0 in range(d0):
+                     for i1 in range(d1):
+                          out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
+            out.resize(sh0[axis_offset:]+sh1[axis_offset:])
+            return out
+     elif isinstance(arg0,escript.Data):
+        if isinstance(arg1,Symbol):
+            return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
+        else:
+            return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
+     else:
+        return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
+ class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
+    """
+    Symbol representing the general tensor product of the transposed of arg0 and arg1
+    """
+    def __init__(self,arg0,arg1,axis_offset=0):
+        """
+        initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
+        @param arg0: first argument
+        @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
+        @param arg1: second argument
+        @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
+        @raise ValueError: inconsistent dimensions of arguments.
+        @note: if both arguments have a spatial dimension, they must equal.
+        """
+        sh_arg0=getShape(arg0)
+        sh_arg1=getShape(arg1)
+        sh01=sh_arg0[:axis_offset]
+        sh10=sh_arg1[:axis_offset]
+        sh0=sh_arg0[axis_offset:]
+        sh1=sh_arg1[axis_offset:]
+        if not sh01==sh10:
+            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)
+        DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
+    def getMyCode(self,argstrs,format="escript"):
+       """
+       returns a program code that can be used to evaluate the symbol.
+       @param argstrs: gives for each argument a string representing the argument for the evaluation.
+       @type argstrs: C{list} of length 2 of C{str}.
+       @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
+       @type format: C{str}
+       @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
+       @rtype: C{str}
+       @raise NotImplementedError: if the requested format is not available
+       """
+       if format=="escript" or format=="str" or format=="text":
+          return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
+       else:
+          raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
+    def substitute(self,argvals):
+       """
+       assigns new values to symbols in the definition of the symbol.
+       The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
+       @param argvals: new values assigned to symbols
+       @type argvals: C{dict} with keywords of type L{Symbol}.
+       @return: result of the substitution process. Operations are executed as much as possible.
+       @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
+       @raise TypeError: if a value for a L{Symbol} cannot be substituted.
+       """
+       if argvals.has_key(self):
+          arg=argvals[self]
+          if self.isAppropriateValue(arg):
+             return arg
+          else:
+             raise TypeError,"%s: new value is not appropriate."%str(self)
+       else:
+          args=self.getSubstitutedArguments(argvals)
+          return generalTransposedTensorProduct(args[0],args[1],args[2])
+ def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
+     "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
+     return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
+ def matrix_transposed_mult(arg0,arg1):
+     """
+     matrix-transposed(matrix) product of the two argument:
+     out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
+     The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
+     The last dimensions of arg0 and arg1 must match.
+     @param arg0: first argument of rank 2
+     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @param arg1: second argument of rank 2
+     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @return: the product of arg0 and the transposed of arg1 at each data point
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     @raise ValueError: if the shapes of the arguments are not appropriate
+     """
+     sh0=getShape(arg0)
+     sh1=getShape(arg1)
+     if not len(sh0)==2 :
+         raise ValueError,"first argument must have rank 2"
+     if not len(sh1)==2 and not len(sh1)==1:
+         raise ValueError,"second argument must have rank 1 or 2"
+     return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
+ def tensor_transposed_mult(arg0,arg1):
+     """
+     the tensor product of the first and the transpose of the second argument
+     for arg0 of rank 2 this is
+     out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
+     and for arg0 of rank 4 this is
+     out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
+     or
+     out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
+     In the first case the the second dimension of arg0 and arg1 must match and
+     in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
+     The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
+     @param arg0: first argument of rank 2 or 4
+     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
+     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
+     @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     """
+     sh0=getShape(arg0)
+     sh1=getShape(arg1)
+     if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
+        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
+     elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
+        return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
+     else:
+         raise ValueError,"first argument must have rank 2 or 4"
+ def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
+     """
+     generalized tensor product of transposed of arg0 and arg1:
+     out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
+     where
+         - s runs through arg0.Shape[:arg0.Rank-axis_offset]
+         - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
+         - t runs through arg1.Shape[arg1.Rank-axis_offset:]
+     The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
+     to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
+     @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 general tensor product of transposed(arg0) and arg1 at each data point.
+     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
+     """
+     if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
+     arg0,arg1=matchType(arg0,arg1)
+     # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
+     if isinstance(arg0,numarray.NumArray):
+        if isinstance(arg1,Symbol):
+            return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
+        else:
+            if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-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)
+            arg0_c=arg0.copy()
+            arg1_c=arg1.copy()
+            sh0,sh1=arg0.shape,arg1.shape
+            d0,d1,d01=1,1,1
+            for i in sh0[:arg0.rank-axis_offset]: d0*=i
+            for i in sh1[:arg1.rank-axis_offset]: d1*=i
+            for i in sh1[arg1.rank-axis_offset:]: d01*=i
+            arg0_c.resize((d0,d01))
+            arg1_c.resize((d1,d01))
+            out=numarray.zeros((d0,d1),numarray.Float64)
+            for i0 in range(d0):
+                     for i1 in range(d1):
+                          out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
+            out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
+            return out
+     elif isinstance(arg0,escript.Data):
+        if isinstance(arg1,Symbol):
+            return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
+        else:
+            return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
+     else:
+        return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
+ class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
+    """
+    Symbol representing the general tensor product of arg0 and the transpose of arg1
+    """
+    def __init__(self,arg0,arg1,axis_offset=0):
+        """
+        initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
+        @param arg0: first argument
+        @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
+        @param arg1: second argument
+        @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
+        @raise ValueError: inconsistent dimensions of arguments.
+        @note: if both arguments have a spatial dimension, they must equal.
+        """
+        sh_arg0=getShape(arg0)
+        sh_arg1=getShape(arg1)
+        sh0=sh_arg0[:len(sh_arg0)-axis_offset]
+        sh01=sh_arg0[len(sh_arg0)-axis_offset:]
+        sh10=sh_arg1[len(sh_arg1)-axis_offset:]
+        sh1=sh_arg1[:len(sh_arg1)-axis_offset]
+        if not sh01==sh10:
+            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)
+        DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
+    def getMyCode(self,argstrs,format="escript"):
+       """
+       returns a program code that can be used to evaluate the symbol.
+       @param argstrs: gives for each argument a string representing the argument for the evaluation.
+       @type argstrs: C{list} of length 2 of C{str}.
+       @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
+       @type format: C{str}
+       @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
+       @rtype: C{str}
+       @raise NotImplementedError: if the requested format is not available
+       """
+       if format=="escript" or format=="str" or format=="text":
+          return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
+       else:
+          raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
+    def substitute(self,argvals):
+       """
+       assigns new values to symbols in the definition of the symbol.
+       The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
+       @param argvals: new values assigned to symbols
+       @type argvals: C{dict} with keywords of type L{Symbol}.
+       @return: result of the substitution process. Operations are executed as much as possible.
+       @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
+       @raise TypeError: if a value for a L{Symbol} cannot be substituted.
+       """
+       if argvals.has_key(self):
+          arg=argvals[self]
+          if self.isAppropriateValue(arg):
+             return arg
+          else:
+             raise TypeError,"%s: new value is not appropriate."%str(self)
+       else:
+          args=self.getSubstitutedArguments(argvals)
+          return generalTensorTransposedProduct(args[0],args[1],args[2])
+ def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
+     "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
+     return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
  #=========================================================
  #  functions dealing with spatial dependency

 Legend:



Removed from v.775
 


changed lines


 
Added in v.1044
 Legend:



Removed from v.775
 


changed lines


 
Added in v.1044
-Removed from v.775
+Added in v.1044

	ViewVC Help
Powered by ViewVC 1.1.26