ufl team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, May 06, 2008 at 10:47:25PM +0200, Martin Sandve Alnæs wrote:
> 2008/5/6 Anders Logg <logg@xxxxxxxxx>:
> > On Tue, May 06, 2008 at 10:38:46AM +0200, Martin Sandve Alnæs wrote:
> >> 2008/5/6 Anders Logg <logg@xxxxxxxxx>:
> >> > On Tue, May 06, 2008 at 10:26:27AM +0200, Martin Sandve Alnæs wrote:
> >> >> 2008/5/6 Anders Logg <logg@xxxxxxxxx>:
> >> >> > On Mon, May 05, 2008 at 10:45:25PM +0200, Martin Sandve Alnæs wrote:
> >> >> >> 2008/5/5 Anders Logg <logg@xxxxxxxxx>:
> >> >> >> > On Mon, May 05, 2008 at 03:34:41PM +0200, Martin Sandve Alnæs wrote:
> >> >> >> >> 2008/5/5 Anders Logg <logg@xxxxxxxxx>:
> >> >> >> >> > On Mon, May 05, 2008 at 03:05:24PM +0200, Martin Sandve Alnæs wrote:
> >> >> >> >> >> 2008/5/5 Anders Logg <logg@xxxxxxxxx>:
> >> >> >> >> >> > On Mon, May 05, 2008 at 02:02:39PM +0200, Martin Sandve Alnæs wrote:
> >> >> >> >> >> >> 2008/5/2 UFL <ufl@xxxxxxxxxx>:
> >> >> >> >> >> >> > 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 <logg@xxxxxxxxx>
> >> >> >> >> >> >> > 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.
> >> >> >> >> >> >>
> >> >> >> >> >> >> I don't agree with this, we want to have e.g. symmetric
> >> >> >> >> >> >> second order tensor elements, which do not fit into this.
> >> >> >> >> >> >
> >> >> >> >> >> > In what way do they not fit?
> >> >> >> >> >> >
> >> >> >> >> >> > 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.
> >> >> >> >> >
> >> >> >> >> > 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.
> >> >> >> >>
> >> >> >> >> I'm not sure, but I feel that this should closely match the resulting
> >> >> >> >> 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.
> >> >> >> >
> >> >> >>
> >> >> >> 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.
> >> >> >
> >> >> >> But if I interpret this correctly as Python list, I don't see the
> >> >> >> advantage of storing this kind of list with None values.
> >> >> >
> >> >> > No, it's nested classes. The list notation was just shorthand (and
> >> >> > confusing).
> >> >> >
> >> >> >> 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:
> >> >> >
> >> >> > How would you store a simple Taylor-Hood element that way?
> >> >>
> >> >> 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).
> >> >
> >> >> > 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.
> >> >>
> >> >> 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.
> >> >
> >> > 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.
> >> >
> >>
> >> 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.
> >
> > Do you want the sub elements of TensorElement and VectorElement to
> > always be scalar?
> 
> In principle, yes, that would feel natural to me.
> (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.

> > 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.
> 
> 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.
> 
> In fact, what you write here is close to this idea:
> 
> >  BDM = FiniteElement("BDM", ...)
> >
> >  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)
> >
> >  print div(u)
> 
> 
> 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.

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

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

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