ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

It could be added with some effort, but I'd prefer to wait until the
redesign of the form language (UFL) is in place.

I'll add tensor-valued functions to the design specification. Any
other suggestions are welcome.

/Anders


On Wed, Sep 27, 2006 at 02:28:39PM +0200, Garth N. Wells wrote:
> Somewhat connected to this discussion, how hard will it be to add
> tensor-valued functions to FFC? This would be a great addition as it is
> tedious to program and a significant source of errors. I can image that
> this is easy from the DOLFIN side since FFC generates the dof mapping.
> 
> Garth
> 
> 
> Anders Logg wrote:
> > 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
> > 
> > _______________________________________________
> > FFC-dev mailing list
> > FFC-dev@xxxxxxxxxx
> > http://www.fenics.org/mailman/listinfo/ffc-dev
> > 
> 
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/ffc-dev