ffc team mailing list archive

[Johan Jansson <jjan@xxxxxxxxxxx>]

On Tue, Sep 26, 2006 at 05:38:17PM +0200, Anders Logg wrote:
> On Tue, Sep 26, 2006 at 05:11:24PM +0200, Dag Lindbo wrote:
> > 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
> 
> 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.
> 

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:

>>> from Numeric import *
>>> x = array([[1, 1, 1]])
>>> matrixmultiply(x, transpose(x))
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