ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Aug 31, 2005 at 07:01:02PM +0200, Johan Jansson wrote:
> 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))

ok, now I see. Yes, this is a tensor (or more correctly a particular
component of a tensor). But it's not indexable since the indices are
already in there (i and j).

You can then declare epsilonv in the same way and do the tensor
contraction by

	    epsilonu * epsilonv

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

Everything in FFC is a scalar, except for vector-valued basis
functions which need to be indexed in the declaration of a form to
produce something that is scalar. More precisely: everything in FFC is
an already indexed tensor of arbitrary rank. So if you do

   u[i].dx(j).dx(k)

then you have a scalar component of a rank 3 tensor, but with an index
(i, j, k) that need not be fixed.

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

When you do tmp = u[i].dx(j), you first pick one component of the
vector-valued basis function u, so then you have a scalar. When you
then take the derivative in the j direction you obtain a new scalar
function.

It's like when you write the same thing on a piece of paper. If you
have a vector x, then x_i is not a vector but component i on the
vector. Sometimes you use the notation (x_i) to denote the vector with
components x_i for i = 1,2,.... Maybe we could add something like this
to FFC, possibly a keyword vec so that

    vec(u[i])

is the thing you want.

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

You can take two vector-valued BasisFunctions u and v and add them to
get a vector-valued Sum:

    u + v

and then you can pick a component:

    (u + v)[i] = u[i] + v[i]

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

No, same explanation as above.

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

They are in a sense matrix-valued already (like u[i].dx(j)) but you
always work with the scalar components.

I hope things are a little more clear now. I need to write this up in
the FFC manual to avoid confusion.

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

/Anders