Subject: [UFL-dev] Evaluating tensor expressions

From: Anders Logg <[email-masked]>
Date: Tue, 21 Oct 2008 17:27:55 +0200

Is there a simple way to evaluate tensor expressions? For example, if

 a = dot(v, u)*dx

is there then a way to convert this to

 v[0]*u[0] + v[1]*u[1] + ...

?

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Tue, 21 Oct 2008 17:34:29 +0200

>  v[0]*u[0] + v[1]*u[1] + ...

I found the following in Dot:

  def as_basic(self, dim, a, b):
      ii = Index()
      aa = a[ii] if (a.rank() == 1) else a[...,ii]
      bb = b[ii] if (b.rank() == 1) else b[ii,...]
      return aa*bb

Why are the arguments dim, a and b required? A Dot already knows its
operands.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 21 Oct 2008 17:41:15 +0200

> Why are the arguments dim, a and b required? A Dot already knows its

Because it is used in an algorithm that first expands its operands,
so a and b here are not necessarily its original operands.
As you see, dim is not used here, but it's used in some other Compounds.

These functions could of course also be placed in the algorithm
in question instead of in each Compound object, like with several
other algorithms. But I think that's a detail we can discuss later,
in context with other similar issues.

-- 
Martin

---

From: Anders Logg <[email-masked]>
Date: Tue, 21 Oct 2008 17:54:46 +0200

> in context with other similar issues.

Yes, if as_basic does not work like a member function, it seems
strange to make it a member function. Maybe it could be made static if
necessary.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 21 Oct 2008 18:03:27 +0200

> strange to make it a member function. Maybe it could be made static if

I'm fine with moving it into the algorithm, I can do that later.

Keeping it (static) in the class simplifies adding new Compound
objects, that's the only reason I can think of now. But this is not
really necessary since they would have to be handled in e.g. the
latex compiler and maybe other places for that to make any sense.

This probably makes no sense to you right now, we should set
up a session where I show how these things fit together.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Tue, 21 Oct 2008 18:07:02 +0200

> up a session where I show how these things fit together.

There are a lot of things that don't make sense to me right now... :-)

So is there any way right now to handle compound expressions, or do I
need to implement it?

Say that I have a = dot(v, u)*dx, how can I expand it into simpler
operations (including only +, *, and [])?

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 21 Oct 2008 18:14:00 +0200

> operations (including only +, *, and [])?

As I just mentioned, the function expand_compounds will expand
dot, inner, grad, div, and several others into index-notation.
Not sure what it is you want beyond that?

Expanding implicit sums and products into explicit sums and
products is not implemented in UFL, since so far I've opted
to do that while converting to swiginac in SFC. It would be nice
to have it in UFL though, but I didn't think you needed that.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Tue, 21 Oct 2008 18:39:00 +0200

> to have it in UFL though, but I didn't think you needed that.

What is the dim argument to expand_compounds, and why is it needed?

Should I look at the operands and find out the appropriate dimension,
and how is that best done?

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Tue, 21 Oct 2008 18:18:51 +0200

> to have it in UFL though, but I didn't think you needed that.

Your explanation of expand_compounds was not addressed to the list and
ended up in another inbox so I missed it... :-) Thanks.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 21 Oct 2008 19:57:01 +0200

> and how is that best done?

As you've discovered, you don't need to use that.

But anyway, dim is the "default dimension", i.e. the spatial dimension.
It's used by as_basic in at least the inverse and cofactor, which
are simply hardcoded with expressions I obtained from ginac...

There are some issues with this dimension / shape thing,
but I'll leave that for another day.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Tue, 21 Oct 2008 20:20:45 +0200

> are simply hardcoded with expressions I obtained from ginac...

Does not sound optimal...

> but I'll leave that for another day.

I'll use dim=None for now and hope we can remove the argument at some
point.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 21 Oct 2008 20:22:48 +0200

> I'll use dim=None for now and hope we can remove the argument at some

You shouldn't need to call it at all, why do you do that?

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Tue, 21 Oct 2008 20:53:43 +0200

> You shouldn't need to call it at all, why do you do that?

What do you mean? I want to call expand_compounds to get from Dot,
Grad etc to index notation. And expand_compounds requires a dim
argument.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 21 Oct 2008 21:11:45 +0200

> Grad etc to index notation. And expand_compounds requires a dim

Sorry, I was thinking that expand_compounds figured out the dim.
Given a Form, you easily know the right dim, but expand_compounds only
takes an expression.

I've moved the as_basic functions now.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Tue, 21 Oct 2008 21:24:22 +0200

> Given a Form, you easily know the right dim, but expand_compounds only

ok, so dim is always just the spatial dimension (as opposed to space
dimension or vector dimension)?

> I've moved the as_basic functions now.

ok!

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 21 Oct 2008 22:01:26 +0200

> dimension or vector dimension)?

What do you mean by the space dimension or vector dimension?
Are you thinking of e.g. 2D problems in 3D space?
I haven't really considered that.

I took a little look at the code again. The reason for dim is to get
the shape of the operand expressions right. The "sinners" are:

In geometry.py:
class FacetNormal(Terminal):
    def shape(self):
        return (DefaultDim,)

Since this is a Terminal, the only way it may know its correct
shape is to have e.g. the polygon as a constructor argument.
But then we can't have a default object "n = FacetNormal()".
As long as this isn't solved, we can't guarantee the invariant
that all shapes are known.

differentiation.py:
class Grad(UFLObject):
    def shape(self):
        return (DefaultDim,) + self._f.shape()

class Curl(UFLObject):
    def shape(self):
        return (DefaultDim,)

These are not terminals, so they might get the spatial
dimension from the subexpression, only... I can't see that
all subexpressions contain a defined spatial dimension.

tensors.py:
class ComponentTensor:
this needs to know the dimensions of the free indices of
its expression, which leads to some unsolved issues of its own,
but this might turn out fine if we can solve everything above.

--
Martin