Subject: [UFL-dev] [HG UFL] Implement recursive initialization of tensor elements.

From: [email-masked] (UFL)
Date: Fri, 2 May 2008 23:16:12 +0200 (CEST)

One or more new changesets pushed to the primary ufl repository.
A short summary of the last three changesets is included below.

changeset:   99:ab1f19b9eb5545605d2651bf95508956aa27ad0b
tag:         tip
user:        Anders Logg <[email-masked]>
date:        Fri May 02 23:16:10 2008 +0200
files:       ufl/basisfunctions.py ufl/finiteelement.py
description:
Implement recursive initialization of tensor elements.
Unit tests pass, but see what you think of it.

changeset:   98:146f65c19d1feca87865470e34878170d7359549
user:        Anders Logg <[email-masked]>
date:        Fri May 02 22:38:27 2008 +0200
files:       ufl/finiteelement.py
description:
Minor fix in pretty printing of mixed elements

changeset:   97:5424a96b780b24a009b823d47d466d64e9d23531
user:        Anders Logg <[email-masked]>
date:        Fri May 02 22:37:13 2008 +0200
files:       ROADMAP ufl/elements.py ufl/finiteelement.py
description:
Change all member data to _foo and add corresponding foo()
functions in finite element classes, most of them collected
in the common FiniteElementBase base class.

----------------------------------------------------------------------
For more details, visit http://www.fenics.org/hg/ufl

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 5 May 2008 14:02:39 +0200

> Unit tests pass, but see what you think of it.

I don't agree with this, we want to have e.g. symmetric
second order tensor elements, which do not fit into this.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 5 May 2008 14:45:19 +0200

> second order tensor elements, which do not fit into this.

In what way do they not fit?

The flag is_symmetric() is still there.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 5 May 2008 15:05:24 +0200

> The flag is_symmetric() is still there.

For a symmetric tensor in 3D, there should be only 6 subelements,
and this recursive initialization will create 9.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 5 May 2008 15:27:58 +0200

> and this recursive initialization will create 9.

Does that matter? Creating an UFL element is very lightweight so it
shouldn't matter that we have a few extra elements. The important
thing is for algorithms (and form compilers) to look at is_symmetric
and use that to optimize.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 5 May 2008 15:34:41 +0200

> and use that to optimize.

I'm not sure, but I feel that this should closely match the resulting
UFC element hierarchy that comes out at the other side.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 5 May 2008 16:27:52 +0200

> UFC element hierarchy that comes out at the other side.

Then what should the hierarchy look like for a symmetric 2x2 tensor
element?

We could do

  [[e, e], [None, e]]

if we like.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 5 May 2008 22:45:25 +0200

> Then what should the hierarchy look like for a symmetric 2x2 tensor

Basically, yes, but I'm not sure what you mean with the []'s?
But if I interpret this correctly as Python list, I don't see the
advantage of storing this kind of list with None values.

I'd rather store a single list of subelements, which can
be handled like in any MixedElement most of the time,
and store the subelement indices like:

2D:
[
 [0, 1],
 [1, 2]
]

3D:
[
 [0, 1, 2],
 [1, 3, 4],
 [2, 4, 5]
]

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Tue, 6 May 2008 10:16:26 +0200

> Basically, yes, but I'm not sure what you mean with the []'s?

I mean that we have a mixed element which has two sub elements,
each consisting of two (scalar) elements.

> advantage of storing this kind of list with None values.

No, it's nested classes. The list notation was just shorthand (and
confusing).

> and store the subelement indices like:

How would you store a simple Taylor-Hood element that way?

I think the natural thing would be to have a MixedElement which has
two sub elements. The first one of these (for u) will be a
VectorElement (which is a special case of a MixedElement) and that
element will have a list of two scalar sub elements. The second one
(for p) will be a FiniteElement (a scalar).

Similarly, one can take a vector-valued BDM-element for u and a scalar
element for p and combine them into a MixedElement. In that case, each
of the two sub elements will be a FiniteElement, but one will be
vector-valued and one scalar.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 6 May 2008 10:26:27 +0200

> How would you store a simple Taylor-Hood element that way?

Irrelevant. A Taylor-Hood element is not a TensorElement.

> vector-valued and one scalar.

I completely agree with that for a general MixedElement.

I just don't agree that a TensorElement should be
represented as a nested hierarchy of MixedElements,
since this will not be simple as intended for symmetric elements.

Actually, I'm still not quite comfortable with VectorElement meaning
"vector of element" instead of "vector valued element", and I don't
think that matches what people would expect. Our different opinions
here may be related to this.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Tue, 6 May 2008 10:31:12 +0200

> Irrelevant. A Taylor-Hood element is not a TensorElement.

No, it's not irrelevant if a TensorElement should be a subclass of
MixedElement (which I think it should since it is a special case).

> here may be related to this.

It's very convenient if one can take any scalar element and use that
to create a vector-valued element, and it's also convenient if one can
have a single base class (MixedElement) which takes care of it all
and have VectorElement and TensorElement as special cases.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 6 May 2008 10:38:46 +0200

> and have VectorElement and TensorElement as special cases.

I never said TensorElement shouldn't be a MixedElement.
I just want it to be a single MixedElement, instead of
containing a nested hierarchy of MixedElements.
The reason is that I don't think symmetric elements will fit
nicely into the hierarchic approach. I suspect using "None"
as you suggested above will lead to all sorts of special cases.

--
Martin

---

From: Kent-Andre Mardal <[email-masked]>
Date: Tue, 06 May 2008 11:13:06 +0200

> as you suggested above will lead to all sorts of special cases.

I think a mixed element is a collection of (typically different)
elements, a vector element is an element with values that are vectors, a
tensor element is an element with values that are tensors. Symmetric
tensors are useful because then appear often.  

You can of course create a 2D vector element as a mixed element
containing two scalar elements. This is not typical, but can be done.
Still this should not make VectorElement a subclass of MixedElement or
vica versa.  

The famous AFW elasticity elements is a mixed element where the first
element is a (possibly symmetric) tensor, the second element is a vector
element and additionally you can have a scalar element if you want. 

Kent

---

From: Anders Logg <[email-masked]>
Date: Tue, 6 May 2008 22:05:01 +0200

> Still this should not make VectorElement a subclass of MixedElement or

I think it's typical. And it's also the way it's implemented in UFC.
At least this is what FFC does; if you create a vector Lagrange
element then you can ask for the sub elements representing each of the
components. This means you can do things like

  w = solve(stokes)

  u = w[0]
  p = w[1]

  ux = u[0]
  uy = u[1]

FIAT used to have vector-valued Lagrange elements but I think they
have been removed since they are not used by FFC anymore. FFC just
takes the scalar elements and adds them up to a mixed element.

-- 
Anders

> _______________________________________________
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev

---

From: Anders Logg <[email-masked]>
Date: Tue, 6 May 2008 22:09:19 +0200

> as you suggested above will lead to all sorts of special cases.

Do you want the sub elements of TensorElement and VectorElement to
always be scalar?

It looks to me like that would prevent us from doing things like
taking 3 BDM elements and combining them into a matrix-valued element
where each row of the matrix is a BDM-function and then operating on
the basis functions like they were matrix-valued functions. I think
Marie does this.

  BDM = FiniteElement("BDM", ...)

  element = TensorElement(BDM, 3)
  
  u = BasisFunction(element)

  print div(u)

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 6 May 2008 22:47:25 +0200

> Do you want the sub elements of TensorElement and VectorElement to

In principle, yes, that would feel natural to me.
(But see below, I'm not too religious about this, really).

> the basis functions like they were matrix-valued functions. I think

No, because we could allow combining multiple vector
valued elements into matrix-valued elements by definition.
But I see that it could make it slightly more difficult to
combine 3 BDM elements and extract each BDM element
as a subelement.

>  element = TensorElement(BDM, 3)

Is this line what you intended? Because here you write TensorElement
in a context I interpret as "tensor/matrix valued element",
which is exactly my point.

>  u = BasisFunction(element)

Personally, I can  _get used to_  VectorElement meaning "vector of elements"
instead of "vector valued element", although in principle I believe
the _intuitive_
meaning will be the opposite for most people. I think we should ask for the
opinion of some non-developers to check this out. Other than that, I won't
object to this anymore.

Also, when it comes to general tensor elements as in "tensors of elements",
your recursive definition is great. So, for non-symmetric tensor elements,
I don't have any problem with them.

Question: Does _anybody_ _really_ care about tensor elements of order >2?
They can always be constructed by manual recursive application
of MixedElement, so I think we can safely ignore them.

In particular, the specification of general tensor elements
of order >2 with particular symmetries becomes pretty involved.
It is probably easier for a user to define his own numbering
for a tensor element of order >2 with symmetries than to care
about any conventions we define (which we won't use).

The user can then take a VectorElement, make a [Basis]Function
on this element, and make a Tensor from components of
this [Basis]Function.

Thus we may reduce our differences to the case of symmetric
tensor elements of order 2 (which is what will be used anyway).
I think we can sort this out without any problems.
This information will be used mostly inside form compilers anyway.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Sat, 10 May 2008 09:00:22 +0200

> (But see below, I'm not too religious about this, really).

I'm starting to think it's a good idea. It would simplify the
implementation.

> opinion of some non-developers to check this out. Other than that, I won't

I don't think there's a big difference, or at least we should try to
minimize the difference. It's a matter of what you can and can not do
with basis functions and elements.

I think you should always be able to pick components of basis
functions (if value rank is nonzero), even for mixed elements:

  element = MixedElement(...)
  v = BasisFunction(element)
  print v[(i, j, k)]

The first index will pick out subelement i, the second will
pick out subelement j of that subelement and then so on until we are
out of subelements, and if any index remains it will pick out
component k of the final (simple) element.

Then sometimes you will be able to pick out subelements and sometimes
not. So for a vector-valued Lagrange element, you will be able to pick
out the subelements, but not for a BDM element since it has no
subelements (it is simple).

> This information will be used mostly inside form compilers anyway.

We might instead want a specific class for this:
SymmetricTensorElement.

Anyway, my suggestion would be that we do the following:

1. MixedElement as it is now but we add __getitem__ which should
return a tuple consisting of element and component associated with the
indexing. For a vector Lagrange element, you would get

  element[(i,)] = (scalar Lagrange element, None) 

For a BDM element, you would get

  element[(i,)] = (BDM element, (,))

Note that this makes it possible to create a matrix-valued element as
a mixed element consisting of three BDM elements and index the basis
functions as matrices.

2. VectorElement and TensorElement only work for scalar elements and
TensorElement accepts an extra argument symmetry which is a mapping
from indices to indices that lays out the symmetric structure. That's
all we need to return the right thing in __getitem__ (which we
overload) and it will support all kinds of strange user-defined
symmetries.

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Sat, 10 May 2008 14:19:09 +0200

> a mixed element consisting of three BDM elements and index the basis

I don't quite understand the component part here. What is it for?
And I was expecting i for the vector Lagrange element.
"All components" can be represented by the built-in
object "...", aka Ellipsis, which means "take all" in a[...].

But maybe it would be clearer to just return the subelement,
and have a function to get the component.

  sub = element[i]
  comp = element.component(i)

We should use element[i], not element[(i,)].

> overload) and it will support all kinds of strange user-defined

Ok.

-- 
Martin

---

From: Anders Logg <[email-masked]>
Date: Sun, 11 May 2008 21:14:42 +0200

>   comp = element.component(i)

What I mean is to provide the one and only information a form compiler
needs from a BasisFunction/Function, namely which simple element it is
associated with and which component (if any) of that element.

I assume here that whenever a BasisFunction/Function appears in a
form, it does so ultimately as a scalar (indexed in some way). At
least this is the case in FFC and I guess in general?

So if a Lagrange function v appears in a form as v[i], we need to know
two things from v[i], namely

  1. It comes from a Lagrange basis function
  2. It is component i of that basis function

Similarly, if a function v from a matrix-valued mixed element
constructed by creating a MixedElement from 3 BDM elments appears
in a form as v[(i, j)], we need to know two things from v[(i, j)],
namely

  1. It comes from a BDM basis function
  2. It is component j of that basis function

So in summary we need to provide something that maps a component
i = (i0, i1, ...) for a BasisFunction/Function of an element (mixed or
simple) to

  1. A simple element
  2. A component of that element

The component will sometimes be equal to the tuple i = (i0, i1, ...)
and sometimes it will be a subset like i' = (i3, i4, ...).

And my suggestion would be to use __getitem__ to return both the
simple element and the tuple.

> We should use element[i], not element[(i,)].

Agree.

-- 
Anders