ffc team mailing list archive

[Johan Jansson <johanjan@xxxxxxxxxxxxxxxx>]

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

...

> 2. Then Johan's suggestion. If you write
> 
>     epsilonu = 0.5 * (u[i].dx(i) * v[j].dx(j))
> 
> then this is a scalar since repeating an Index twice means summation
> over that Index. I think the tensor notation has proven to be quite
> powerful so I would like to keep it.

I'm very sorry, but my expression was totally wrong. I copied it from
the wrong place, and I didn't catch it because it looked pretty
similar. What I meant was:

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

Here the indices are not repeated, so it should not be a scalar. But
it's still not indexable in FFC. Is it meant to be indexable? Are
Products and Sums also supposed to be tensors?

Here's a minimal example:

# Tensor test

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

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

gradu = u[i].dx(j)

tmp = gradu[0]

a = tmp * v[0] * dx

When I compile this, I get:

  File "Foo.py", line 36, in ?
    tmp = gradu[0]
  File "/usr/lib/python2.4/site-packages/ffc/compiler/algebra.py", line 383, in __getitem__
    w.basisfunctions[0] = w.basisfunctions[0][component]
  File "/usr/lib/python2.4/site-packages/ffc/compiler/algebra.py", line 209, in __getitem__
    raise RuntimeError, "Cannot pick component of scalar BasisFunction."
RuntimeError: Cannot pick component of scalar BasisFunction.

Perhaps it's just a bug in FFC? I haven't been able to follow the
error clearly. I see now that Product and Sum do have an indexing
operator, so perhaps I was just confused about this. If they are
indeed tensors, then we can just as well define tensor operators for
them.

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

Here you construct a nested list structure which is a sort of tensor
representation though (like Numeric.array). Shouldn't this work (i and
j are Index)?:

def epsilon(u):
    return u[i].dx(j) + u[j].dx(i)

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

It would also be good to have matrix operators, but that requires
matrix-valued algebra objects. Matrix multiplication, gradient of a
vector etc. But as I mentioned, as soon as indexing works, then the
rest of the operators should be trivial to implement.

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

Yes, I agree with this.

  Johan