Subject: [FFC-dev] Tensors

From: Johan Jansson <[email-masked]>
Date: Tue, 23 Aug 2005 14:19:28 +0200

Hi,

I'm at a point where tensor (matrix) notation would be very useful in
FFC. Has anyone given any thought to what needs to be done? I can work
on implementing it, I need some exercise in Python and to learn more
about FFC. It would be good to discuss which features are needed
here on this list.

Right now I'm thinking that all that should be necessary is a matrix
indexing operator: [a, b] (perhaps general tensor indexing can wait,
but presumably it would extend naturally: [a, b, c] etc.), which would
be very similar to the vector indexing operator: [a]. There also needs
to be a way of specifying the dimensions for the components of the
tensor. Perhaps slicing would also be useful, i.e. the [:, b] and [a,
:] operators.

The actual tensor operators (matrix multiplication etc.) could then be
user-defined, at least at this stage.

What do you think? The rest of FFC and FIAT shouldn't have to know
about tensors, an M x N tensor-valued element should just look like an
M * N vector-valued element to them. I haven't looked much at the code
yet, but hopefully it should be straightforward to do it this way.

  Johan

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 24 Aug 2005 13:32:19 +0200

Hi Johan,

For me, tensor notation would be really useful. 

Related to the issue of taking slices and indexing, I had some exchanges with
Anders a while  ago about the possibility of being able to specify the range of
the summation when using the index summation convention. This would be a big
help when examining coupled problems. At the moment I'm not able (or don't know
how) to use the index notation for these problems, so I'm writing out the
bilinear form for explicitly. As an example, I'm  looking at a coupled
elasticity-diffusion equation at the moment. In 2D (3 dof per node), part of th
bilinear form looks like

a = ( lmbda * (u[0].dx(0) + u[1].dx(1))*(v[0].dx(0) + v[1].dx(1))  + \
       mu * ( \
        (u[0].dx(0)+u[0].dx(0))*(v[0].dx(0)) + \
        (u[0].dx(1)+u[1].dx(0))*(v[0].dx(1)) + \
        (u[1].dx(0)+u[0].dx(1))*(v[1].dx(0)) + \
        (u[1].dx(1)+u[1].dx(1))*(v[1].dx(1)) ) )* dx + \
     -Beta*( u[2]*v[0].dx(0) + u[2]*v[1].dx(1) )* dx

For 3D, I need to add a bunch of terms. What I would like to able to do is use
something like:

a = ( lmbda * u[i:0,n].dx(i:0,n) * v[j:0,n].dx(j:0,n) + \
  mu * ( u[i:0,n].dx(j:0,n) + u[j:0,n].dx(i:0,n) ) * v[i:0,n].dx(j:0,n) )* dx + \
 -Beta * u[n+1]*v[i:0,n].dx(i:0,n) * dx
 
where n is the range of the summation. Setting n = 1 would then be for 2D, and n
 = 2 for 3D. It might then be an idea to make the spatial dimension available 
within ffc.

Garth

Quoting Johan Jansson <[email-masked]>:

> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev
>

---

From: Johan Jansson <[email-masked]>
Date: Mon, 29 Aug 2005 15:56:56 +0200

> For me, tensor notation would be really useful. 

Hi!

I've experimented a bit with different approaches to tensor notation,
and I've found one that works pretty well, some feedback would be
nice. The approach I've gotten to work is to simply define matrices of
objects in the FFC algebra, I've done it using multiarrays from the
Numeric module, but any capable multiarray/matrix class could probably
be used, or we could write our own if we need special capabilities.

What I would like to be able to do is write:

epsilonu = 0.5 * (u[i].dx(i) * v[j].dx(j))

and then treat epsilonu as a tensor. But as far as I can see, epsilonu
is not a tensor in the FFC algebra. At least indexing is not yet
supported for such objects, so I cannot write epsilonu[1, 2] or even
epsilonu[1]. I'm not sure epsilonu is supposed to fit the tensor
concept here.

What I did instead is to explicitly construct a tensor (matrix):

from Numeric import *
from math import *

def epsilon(u):
    eu = ffczeros(3, 3)

    for ii in range(3):
        for jj in range(3):
            eu[ii, jj] = 0.5 * (u[ii].dx(jj) + u[jj].dx(ii))

    return eu

also defining the function ffczeros():

global zero
zero = 0.0 * Constant()

def ffczeros(m, n):
    A = array((zero))
    A = resize(A, (m, n))
    A *= 0.0

    return A

array() and resize() are Numeric.array() and Numeric.resize(). The
array in Numeric is not really made for arbitrary types, it seems to
only officially support builtin types, but this works well
enough. There is an array class in numarray which does support general
types, but it is not well developed. As I mentioned, we can write our
own matrix/tensor class if needed.

I can now write:

epsilonu = epsilon(u)

to achieve what I want.

It's then very easy to write forms using the same notation as in the
original PDEs. For example, here is the E() linear transformation from
linear elasticity:

def E(e, lmbda, mu):
    Ee = 2.0 * multiply(mu, e)
    for ii in range(3):
        for jj in range(3):
            if ii == jj:
                Ee[ii, jj] += lmbda * trace(e)

    return Ee

I would have liked to write it like so:

def E(e, lmbda, mu):
    Ee = 2.0 * multiply(mu, e) + lmbda * trace(e) * ffcidentity(3)
    return Ee

but I have not yet succeeded in constructing an identity matrix of FFC
algebra objects. It would be nice to have a zero and identity element
in the FFC algebra, they should just be a Constant() with a 0.0 or 1.0
value. I'm not really clear on the Python operators yet either, there
seems to be some ambiguity between when an operator is used from FFC
or Numeric.array.

I have attached the elasticity form from the DOLFIN module written in
the new notation (Elasticity-alternative.form), I've verified that it
produces the same result, so it works. FFC seems to have some problems
with namespaces though, so you need to compile it like so:

ffc Elasticity-alternative.form
python Elasticity-alternative.py

> Related to the issue of taking slices and indexing, I had some
> exchanges with Anders a while ago about the possibility of being
> able to specify the range of the summation when using the index
> summation convention. This would be a big help when examining
> coupled problems. At the moment I'm not able (or don't know how) to
> use the index notation for these problems, so I'm writing out the
> bilinear form for explicitly. As an example, I'm looking at a
> coupled elasticity-diffusion equation at the moment. In 2D (3 dof
> per node), part of th bilinear form looks like

> For 3D, I need to add a bunch of terms. What I would like to able to
> do is use something like:

> a = ( lmbda * u[i:0,n].dx(i:0,n) * v[j:0,n].dx(j:0,n) + mu * (
> u[i:0,n].dx(j:0,n) + u[j:0,n].dx(i:0,n) ) * v[i:0,n].dx(j:0,n) )* dx
> + -Beta * u[n+1]*v[i:0,n].dx(i:0,n) * dx

> where n is the range of the summation. Setting n = 1 would then be
> for 2D, and n = 2 for 3D. It might then be an idea to make the
> spatial dimension available within ffc.  Garth

I think this could be done much simpler with the separate
tensor/matrix class approach. You would do:

def epsilon(u, n):
    eu = ffczeros(n, n)

    for ii in range(n):
        for jj in range(n):
            eu[ii, jj] = 0.5 * (u[ii].dx(jj) + u[jj].dx(ii))

    return eu

and then:

a = epsilon(u, n) * ... -Beta * u[n+1] * ...

Since slicing is supported by Numeric.array, and would be supported by
any other tensor/matrix representation we would choose, it would not
be necessary to add it to the FFC algebra. It would make the matrix
initialization much more compact though, if we could just write:

epsilonu = matrix(u[i:0,n].dx(j:0,n) + u[i:0,n].dx(j:0,n))

(or some other slicing convention), but perhaps it's best to keep the
FFC algebra as simple as possible.

Please tell me what you think about this, if it seems to cover the
necessary functionality etc. I'm going to try to use it for some other
problems and see if it's usable.

  Johan