ffc team mailing list archive

[Johan Jansson <johanjan@xxxxxxxxxxxxxxxx>]

On Wed, Aug 31, 2005 at 10:35:04AM -0500, Anders Logg wrote:

...

> On the other hand, you can do everything you need by indexing basis
> functions and derivatives by integers (not Index) as you have noted.
> This allows a user to define arbitrary operators as you have done:
> 
> def epsilon(u):
>     return [ [u[i].dx(j) + u[j].dx(i) for i in range(d)] for j in range(d) ]
> 
> The only thing that needs to be added to FFC is then the identity as
> far as I can see.
> 
> I'm planning to add a list of predefined operators and operands to FFC
> to simplify and beautify the notation anyway, so this would fit well
> into that category.
> 
> Operators that will be added include grad(), rot() and div().
> 
> Operands that will be added include h (mesh size) and d (space dimension).
> 
> Any other suggestions?
> 
> All these can be added on top of the existing machinery (I think). As
> far as possible, I would like to keep the basic layer simple and add
> complexity by combining the simple primitives.
> 
> /Anders
> 

I think the tensor notation is converging to a natural choice of
representation. I've tested the Numeric approach a bit more, and it's
not actually necessary to have a zero and identity member in the FFC
algebra, the values 0.0 and 1.0 seem to be fine (so it's possible to
use the Numeric functions zeros() and identity() to construct zero and
identity matrices). Some languages have a dynamic type system, but it
seems Python goes even further than that, which is good for us in this
case at least.

I can now write the elasticity form like this (element definitions
etc. not included):

def grad(u):
    return [ [u[i].dx(j) for i in range(d)] for j in range(d) ]

def epsilon(u):
    return 0.5 * (grad(u) + transpose(grad(u)))

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

    return Ee

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

a = vdot(sigma, epsilon(v)) * dx

All that's left for the notation to be identical (as identical as
ASCII can be) to written form is to move to operators instead of
functions for multiply(), vdot() etc. I'm not convinced that's a good
idea though, I prefer this notation, since it's unambigous. A "*"
operator can mean any kind of product. Since there are at least two
products which make sense for matrices in this context (and probably
more in other contexts), I don't think we have enough symbols. It
always takes me a few minutes to figure out what the dot operator
means in an article, perhaps I have finally found a remedy in this
notation convention. It will lead to an excessive amount of
parantheses though, but functional language programmers seem to be
able to handle that just fine (i.e. Lisp).

It also doesn't seem necessary to explicitly construct a Numeric.array
object using the array() function, the list representation works
fine. It does seem necessary for the Numeric.array operators to work
though, so perhaps it would be good to define grad() as:

def grad(u):
    return array(( [ [u[i].dx(j) for i in range(d)] for j in range(d) ] ))

Then we can use the slicing operators (and other operators) from
Numeric.array. For example, to extract the first 2 columns of grad(u),
we can write:

grad(u)[:, 0:2]

Perhaps this can be generalized even further by simply defining a
vec() function Anders mentioned before, and extending the concept to
matrix() and tensor(), which would simply construct a Numeric.array
from an FFC algebra object, i.e.:

i = Index()
j = Index()

gradu = matrix(u[i].dx(j))

which would be a much cleaner fit with matrix-valued elements. I.e. a
matrix-valued function "stress" could be manipulated in the same way:

i = Index()
j = Index()

sigma = matrix(stress[i, j])

This would mean fewer special cases, which I think is almost always
good. I still haven't grasped the Index concept fully though, so
perhaps there are other better ways to accomplish this.

If we need tensors of rank > 2, then we would have to replace
Numeric.array with something else. Is this relevant? Is there anyone
who would actually need that functionality?

  Johan