ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Sep 26, 2006 at 05:49:31PM +0200, Johan Jansson wrote:
> 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
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> 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