ufl team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

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