ufl team mailing list archive

[Kent-Andre Mardal <kent-and@xxxxxxxxx>]

Difficult subject, in particular when it comes to matrices. 
I think we should let * and dot be the same, I mean dot is just a word
for *, right ? 

The meaning of *should be the same as in Matlab, i.e. 
A*B is not the contraction, but the matrix-matrix product. 

Kent


man, 07.04.2008 kl. 11.51 +0200, skrev Martin Sandve Alnæs:
> I'd like to hear opinions of how * should behave in different contexts.
> 
> To begin with, we will of course agree on the meaning of * between
> scalar expressions, which may or may not have free indices. Repeated
> free indices will be interpreted as implicit summation. Single free
> indices are combined to define a larger set of free indices. This is
> already implemented. (Some of the implicit summation is not, but I
> have it on paper).
> 
> We may perhaps also agree that multiplying a matrix (rank 2 tensor)
> with a vector (rank 1 tensor) means contraction over the last index of
> the matrix and the index of the vector, which is equivalent with the
> dot product. This will allow f.ex. Poissons equation with a conduction
> tensor (M Du : Dv) to be written like inner(M*grad(u), grad(v)). Ok?
> 
> But if so, is this to be a special case, or do we simply define * to be dot?
> This will eventually lead to (let u,v be vectors, A,B,C matrices, a scalar):
> 
> Very natural notation:
> u = A*v = dot(A, v) = "A v" in linear algebra notation
> C = A*B = dot(A, B) = "A B" in linear algebra notation
> 
> Feels natural to me:
> u = v*A = dot(v, A) = "v^T A" in linear algebra notation
> a = u*A*v = (u*A)*v = u*(A*v)
>   = dot(u, A*v) = dot(u*A, v)
>   = dot(dot(u, A), v) = dot(u, dot(A, v)) = "u^T A v" in linear algebra notation
> 
> However, keeping * and dot separate may make things simpler in other
> ways in the implementation. And easy code understanding isn't always
> the same as shorter expressions. So what do you think?
> 
> Note that dot(f, g) is defined differently in UFL than in FFC:
>     dot(f, g) == sum(f[...,k] * g[k, ...] for k in range(dim))
> where ... means the obvious for a tensor of rank > 1.
> Using ... to "skip indices" is also implemented.
> The contraction A:B is instead implemented as inner(A, B).
> 
> Another example, Navier-Stokes with explicit convection term (freely
> from memory):
> a = u - w + dt * ( w * grad(u) * v - nu*inner(grad(u), grad(v)) +
> p*div(v) - f*v )
> vs.
> a = u - w + dt * ( dot(dot(w, grad(u)), v) - nu*inner(grad(u),
> grad(v)) + p*div(v) - dot(f, v) )
> 
> --
> Martin
> _______________________________________________
> UFL-dev mailing list
> UFL-dev@xxxxxxxxxx
> http://fenics.org/mailman/listinfo/ufl-dev