ufl team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

2008/5/10 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.

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

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

Ok.

-- 
Martin