Subject: [FFC-dev] outer products

From: "Dag Lindbo" <[email-masked]>
Date: Tue, 26 Sep 2006 17:11:24 +0200 (MEST)

Hi all,

As I've been told, it shoud be possible in FFC to compute an outer product
between a vector valued function and its transpose:

T = mult(vec(n), transp(vec(n)))

a = (...)
L = dot(v,div(T))

Here, n is on a 'Vector Lagrange' element and so is v (test func).

I'm a bit concerned that it's not working: If I put the transpose on the
first term instead I should get something totally different (a rank 0
tensor instead of rank 2, I guess). But when I compute in DOLFIN,
impelmenting the two variants I get the same results.

Have I toatally misunderstood something here?

Regards,
Dag Lindbo

---

From: Anders Logg <[email-masked]>
Date: Tue, 26 Sep 2006 17:38:17 +0200

> Have I toatally misunderstood something here?

Yes, this seems to be a bug. Transpose and matrix multiplications are
handled in FFC by the Python package Numeric and it does not look like
Numeric is aware of the transpose of a vector.

You can try this out yourself in Python:

>>> from Numeric import *
>>> x = array([1,1,1])
>>> matrixmultiply(transpose(x), x)
3
>>> matrixmultiply(x, transpose(x))
3

In both cases, the scalar product is computed and not the outer
product.

One option would be to add an operator outer() in FFC. If you like,
take a look in the file operators.py and see if you can get something
that works. You can either use FFC directly from within Python:

from ffc import *

or put things like this in a .form file:

element = FiniteElement("Vector Lagrange", "triangle", 1)
v = BasisFunction(element)
print mult(transp(vec(v)), vec(v))

and then compile with FFC as usual.

The set of operators defined in operators.py is not as robust and
complete as I would hope and my plan is to replace the form language
with something better (UFL). In the meantime, I'll be happy to accept
any fixes to operators.py. 

/Anders

---

From: Johan Jansson <[email-masked]>
Date: Tue, 26 Sep 2006 17:49:31 +0200

> In both cases, the scalar product is computed and not the outer

I think the problem is that vec() produces a rank-1 tensor, while
transpose() is not defined for rank-1 tensors (it acts like an
identity operator I guess).

What you want is a matrix() operator which produces a rank-2
tensor. Then you can do:

>>> x = array([[1, 1, 1]])
array([       [3]])
>>> matrixmultiply(transpose(x), x)
array([[1, 1, 1],
       [1, 1, 1],
       [1, 1, 1]])

I've implemented this manually in the elasticity-updated forms in
DOLFIN, see for example:

src/modules/elasticity-updated/dolfin/ElasticityUpdated.form

  Johan

---

From: Anders Logg <[email-masked]>
Date: Tue, 26 Sep 2006 17:53:51 +0200

On Tue, Sep 26, 2006 at 05:49:31PM +0200, Johan Jansson wrote:
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev

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

---

From: Johan Jansson <[email-masked]>
Date: Tue, 26 Sep 2006 18:00:46 +0200

On Tue, Sep 26, 2006 at 05:53:51PM +0200, Anders Logg wrote:

...

> make vec() return this directly so it works like a column vector?

The second alternative is probably best. I think mat() should be
reserved for matrix-valued functions, to perform the equivalent of
vec().

  Johan

---

From: Anders Logg <[email-masked]>
Date: Tue, 26 Sep 2006 19:21:22 +0200

> reserved for matrix-valued functions, to perform the equivalent of

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.

/Anders

---

From: "Dag Lindbo" <[email-masked]>
Date: Tue, 26 Sep 2006 19:33:27 +0200 (MEST)

> 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

---

From: Anders Logg <[email-masked]>
Date: Tue, 26 Sep 2006 19:37:00 +0200

> might (if I'm sucessful) be enough until the extension of the language is

ok.

/Anders

---

From: "Dag Lindbo" <[email-masked]>
Date: Wed, 27 Sep 2006 09:45:49 +0200 (MEST)

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:
(...)
f1      = Function(element)
f2      = Function(element)

outer(f1,f2) returns the outer product f1'*f2, a rank-2 tensor.

//Dag Lindbo

---

From: Anders Logg <[email-masked]>
Date: Wed, 27 Sep 2006 09:49:22 +0200

Looks good! I'll add it later today.

/Anders

> """This module extends the form algebra with a collection of operators
> based on the basic form algebra operations."""
> 
> __author__ = "Anders Logg ([email-masked])"
> __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)))

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 27 Sep 2006 14:28:39 +0200

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

---

From: Anders Logg <[email-masked]>
Date: Wed, 27 Sep 2006 19:21:50 +0200

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