ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Looks good! I'll add it later today.

/Anders


On Wed, Sep 27, 2006 at 09:45:49AM +0200, Dag Lindbo wrote:
> Hi,
> 
> I've written a simple operator for the outer prduct as disgussed
> yesterday. The description of how to create a patch in the manual is not
> producing a sensible patch ('make clean' does not do it for python).
> 
> So, I attach my version of 'operators.py'.
> 
> Description:
> element = FiniteElement("Vector Lagrange", "triangle", 1)
> (...)
> f1      = Function(element)
> f2      = Function(element)
> 
> outer(f1,f2) returns the outer product f1'*f2, a rank-2 tensor.
> 
> //Dag Lindbo
> 
> > On Tue, Sep 26, 2006 at 07:33:27PM +0200, Dag Lindbo wrote:
> >> > On Tue, Sep 26, 2006 at 06:00:46PM +0200, Johan Jansson wrote:
> >> >> On Tue, Sep 26, 2006 at 05:53:51PM +0200, Anders Logg wrote:
> >> >>
> >> >> ...
> >> >>
> >> >> > Good point. The following example seems to work:
> >> >> >
> >> >> >   element = FiniteElement("Vector Lagrange", "triangle", 1)
> >> >> >   v = BasisFunction(element)
> >> >> >   print mult(transp([vec(v)]), [vec(v)])
> >> >> >   print mult([vec(v)], transp([vec(v)]))
> >> >> >
> >> >> > Should we add a new operator mat() that returns [vec()] or should
> >> we
> >> >> > make vec() return this directly so it works like a column vector?
> >> >> >
> >> >> >
> >> >> > /Anders
> >> >>
> >> >> The second alternative is probably best. I think mat() should be
> >> >> reserved for matrix-valued functions, to perform the equivalent of
> >> >> vec().
> >> >>
> >> >>   Johan
> >> >
> >> > Sounds good, but this will require some work. A number of other
> >> > operators in operators.py will need to be updated correspondingly.
> >> >
> >> > If anyone is willing to try, you're more than welcome. If not, I'll
> >> > wait until the reimplementation (and extension) of the form language.
> >>
> >> Tomrrow morning I'll write a simple operator 'outer(vec(n),vec(n))' that
> >> might (if I'm sucessful) be enough until the extension of the language
> >> is
> >> complete.
> >>
> >> /Dag
> >
> > ok.
> >
> > /Anders
> >

> """This module extends the form algebra with a collection of operators
> based on the basic form algebra operations."""
> 
> __author__ = "Anders Logg (logg@xxxxxxxxx)"
> __date__ = "2005-09-07 -- 2005-12-20"
> __copyright__ = "Copyright (C) 2005-2006 Anders Logg"
> __license__  = "GNU GPL Version 2"
> 
> # Modified by Ola Skavhaug, 2005
> 
> # Python modules
> import sys
> import Numeric
> 
> # FFC common modules
> sys.path.append("../../")
> from ffc.common.exceptions import *
> 
> # FFC compiler modules
> from index import *
> from algebra import *
> from projection import *
> from finiteelement import *
> 
> def Identity(n):
>     "Return identity matrix of given size."
>     # Let Numeric handle the identity
>     return Numeric.identity(n)
> 
> def rank(v):
>     "Return rank for given object."
>     if isinstance(v, BasisFunction):
>         return v.element.rank() - len(v.component)
>     elif isinstance(v, Product):
>         return rank(v.basisfunctions[0])
>     elif isinstance(v, Sum):
>         return rank(v.products[0])
>     elif isinstance(v, Function):
>         return rank(Sum(v))
>     else:
>         return Numeric.rank(v)
>     return 0
> 
> def vec(v):
>     "Create vector of scalar functions from given vector-valued function."
>     # Check if we already have a vector
>     if isinstance(v, list):
>         return v
>     # Check if we have an element of the algebra
>     if isinstance(v, Element):
>         # Check that we have a vector
>         if not rank(v) == 1:
>             raise FormError, (v, "Cannot create vector from scalar expression.")
>         # Get vector dimension
>         n = __tensordim(v, 0)
>         # Create list of scalar components
>         return [v[i] for i in range(n)]        
>     # Let Numeric handle the conversion
>     if isinstance(v, Numeric.ArrayType) and len(v.shape) == 1:
>         return v.tolist()
>     # Unable to find a proper conversion
>     raise FormError, (v, "Unable to convert given expression to a vector,")
> 
> def dot(v, w):
>     "Return scalar product of given functions."
>     # Check ranks
>     if rank(v) == rank(w) == 1:
>         # Check dimensions
>         if not len(v) == len(w):
>             raise FormError, ((v, w), "Dimensions don't match for scalar product.")
>         # Use index notation if possible
>         if isinstance(v, Element) and isinstance(w, Element):
>             i = Index()
>             return v[i]*w[i]
>         # Otherwise, use Numeric.dot
>         return Numeric.dot(vec(v), vec(w))
>     elif rank(v) == rank(w) == 2:
>         # Check dimensions
>         if not len(v) == len(w):
>             raise FormError, ((v, w), "Dimensions don't match for scalar product.")
>         # Compute dot product (:) of matrices
>         return Numeric.sum([v[i][j]*w[i][j] for i in range(len(v)) for j in range(len(v[i]))])
> 
> def cross(v, w):
>     "Return cross product of given functions."
>     # Check dimensions
>     if not len(v) == len(w):
>         raise FormError, ((v, w), "Cross product only defined for vectors in R^3.")
>     # Compute cross product
>     return [v[1]*w[2] - v[2]*w[1], v[2]*w[0] - v[0]*w[2], v[0]*w[1] - v[1]*w[0]]
> 
> def trace(v):
>     "Return trace of given matrix"
>     # Let Numeric handle the trace
>     return Numeric.trace(v)
> 
> def transp(v):
>     "Return transpose of given matrix."
>     # Let Numeric handle the transpose."
>     return Numeric.transpose(v)
> 
> def mult(v, w):
>     "Compute matrix-matrix product of given matrices."
>     # First, convert to Numeric.array (safe for both array and list arguments)
>     vv = Numeric.array(v)
>     ww = Numeric.array(w)
>     if len(vv.shape) == 0 or len(ww.shape) == 0:
>         # One argument is a scalar
>         return vv*ww
>     if len(vv.shape) == len(ww.shape) == 1:
>         # Vector times vector
>         return Numeric.multiply(vv, ww) 
>     elif len(vv.shape) == 2 and (len(ww.shape) == 1 or len(ww.shape) == 2):
>         # Matvec or matmat product, use matrixmultiply instead
>         return Numeric.matrixmultiply(vv, ww)
>     else:
>         raise FormError, ((v, w), "Dimensions don't match for multiplication.")
> 
> def D(v, i):
>     "Return derivative of v in given coordinate direction."
>     # Use member function dx() if possible
>     if isinstance(v, Element):
>         return v.dx(i)
>     # Otherwise, apply to each component
>     return [D(v[j], i) for j in range(len(v))]
>     
> def grad(v):
>     "Return gradient of given function."
>     # Get shape dimension
>     d = __shapedim(v)
>     # Check if we have a vector
>     if rank(v) == 1:
>         return [ [D(v[i], j) for j in range(d)] for i in range(len(v)) ]
>     # Otherwise assume we have a scalar
>     return [D(v, i) for i in range(d)]
> 
> def div(v):
>     "Return divergence of given function."
>     # Use index notation if possible
>     if isinstance(v, Element):
>         i = Index()
>         return v[i].dx(i)
>     # Otherwise, use Numeric.sum
>     return Numeric.sum([D(v[i], i) for i in range(len(v))])
> 
> def rot(v):
>     "Return rotation of given function."
>     # Check dimensions
>     if not len(v) == __shapedim(v) == 3:
>         raise FormError, (v, "Rotation only defined for v : R^3 --> R^3")
>     # Compute rotation
>     return [D(v[2], 1) - D(v[1], 2), D(v[0], 2) - D(v[2], 0), D(v[1], 0) - D(v[0], 1)]
> 
> def curl(v):
>     "Alternative name for rot."
>     return rot(v)
> 
> def mean(v):
>     "Return mean value of given Function (projection onto piecewise constants)."
>     # Check that we got a Function
>     if not isinstance(v, Function):
>         raise FormError, (v, "Mean values are only supported for Functions.")
>     # Different projections needed for scalar and vector-valued elements
>     element = v.e0
>     if element.rank() == 0:
>         P0 = FiniteElement("Discontinuous Lagrange", element.shape_str, 0)
>         pi = Projection(P0)
>         return pi(v)
>     else:
>         P0 = FiniteElement("Discontinuous vector Lagrange", element.shape_str, 0, element.tensordim(0))
>         pi = Projection(P0)
>         return pi(v)
> 
> def outer(v,w):
>     "Return outer product of vector valued functions, p = v'*w"
>     # Check that we got a Function
>     if not isinstance(v, Function):
>         raise FormError, (v, "Outer products are only defined for Functions.")
>     if not isinstance(w, Function):
>         raise FormError, (w, "Outer products are only defined for Functions.")
>     if not len(v) == len(w):
>         raise FormError, ((v, w),"Invalid operand dims in outer product")
>     
>     vv = vec(v)
>     ww = vec(w)
>     
>     return mult(transp([vv]),[ww])
>     
> def __shapedim(v):
>     "Return shape dimension for given object."
>     if isinstance(v, list):
>         # Check that all components have the same shape dimension
>         for i in range(len(v) - 1):
>             if not __shapedim(v[i]) == __shapedim(v[i + 1]):
>                 raise FormError, (v, "Components have different shape dimensions.")
>         # Return length of first term
>         return __shapedim(v[0])
>     elif isinstance(v, BasisFunction):
>         return v.element.shapedim()
>     elif isinstance(v, Product):
>         return __shapedim(v.basisfunctions[0])
>     elif isinstance(v, Sum):
>         return __shapedim(v.products[0])
>     elif isinstance(v, Function):
>         return __shapedim(Sum(v))
>     else:
>         raise FormError, (v, "Shape dimension is not defined for given expression.")
>     return 0
> 
> def __tensordim(v, i):
>     "Return size of given dimension for given object."
>     if i < 0 or i >= rank(v):
>         raise FormError, ((v, i), "Tensor dimension out of range.")
>     if isinstance(v, BasisFunction):
>         return v.element.tensordim(i + len(v.component))
>     elif isinstance(v, Product):
>         return __tensordim(v.basisfunctions[0], i)
>     elif isinstance(v, Sum):
>         return __tensordim(v.products[0], i)
>     elif isinstance(v, Function):
>         return __tensordim(Sum(v), i)
>     else:
>         raise FormError, ((v, i), "Tensor dimension is not defined for given expression.")
>     return 0
> 
> if __name__ == "__main__":
> 
>     scalar = FiniteElement("Lagrange", "tetrahedron", 2)
>     vector = FiniteElement("Vector Lagrange", "tetrahedron", 2)
> 
>     i = Index()
> 
>     v = BasisFunction(scalar)
>     u = BasisFunction(scalar)
>     w = Function(scalar)
> 
>     V = BasisFunction(vector)
>     U = BasisFunction(vector)
>     W = Function(vector)
>     
>     i = Index()
>     j = Index()
> 
>     dx = Integral()
> 
>     print dot(grad(v), grad(u))*dx
>     print vec(U)
>     print dot(U, V)
>     print dot(vec(V), vec(U))
>     print dot(U, grad(v))
>     print div(U)
>     print dot(rot(V), rot(U))
>     print div(grad(dot(rot(V), U)))*dx
>     print cross(V, U)
>     print trace(mult(Identity(len(V)), grad(V)))
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/ffc-dev