ufl team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Sat, May 10, 2008 at 02:19:09PM +0200, Martin Sandve Alnæs wrote:
> 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)

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


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