ffc team mailing list archive

[Johan Jansson <johanjan@xxxxxxxxxxxxxxxx>]

On Wed, Aug 24, 2005 at 01:32:19PM +0200, Garth N. Wells wrote:
> Hi Johan,
> 
> 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

> 
> 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

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

# Copyright (c) 2004 Anders Logg (logg@xxxxxxxxx)
# Licensed under the GNU GPL Version 2
#
# The bilinear form for a mass matrix.
#
# Compile this form with FFC: ffc Mass.form

from Numeric import *
from math import *

element = FiniteElement("Vector Lagrange", "tetrahedron", 1)

v = BasisFunction(element)
u = BasisFunction(element)

lmbda = Constant() # Lame coefficient
mu    = Constant() # Lame coefficient

f = Function(element) # Source

global zero
zero = 0.0 * Constant()

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

def ffcidentity(m):
    A = myzeros(m, m)

    for ii in range(m):
        for jj in range(m):
            if ii == jj:
                A[ii, jj] = 1.0
                    
    return A

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


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

sigma = E(epsilon(u), lmbda, mu)

a = vdot(sigma, epsilon(v)) * dx
L = f[i] * v[i] * dx