Distinguish enriched and union element in UFL

[Jan Blechta]

FIAT pull request #27 provided a new implementation of enriched element which is nodal. (There is still an fall-back option of using the old non-nodal implementaion in the case that provided subelements are not nodal.) It shows up that this has serious consequences for form compilation. Two principally different options are:

1. Nodal enriched element. Nodes are concatenated, basis functions are reorthogonalized for nodality. Hence, the result is a new element without any structure which could be exploited within form compilation.

2. Enriched element without nodes. Basis functions are concatenated, nodes are dropped. Structure (mixed or enriched) of subelements is exploited during form compilation because basis is hierarchically nested.

This should be reflected in the form language having:

1. EnrichedElement meaning nodal enriched element,

2. UnionElement which would mean the current, nodeless implementation.

[Lawrence Mitchell]

This sounds like a bad idea to me. UFL knows nothing about basis functions, and it would be good to keep it that way. This smells like exposing internal implementation details in FIAT to UFL, which the user must then know. When do I use EnrichedElement over UnionElement and vice versa? Well, I have to go and look in the FIAT source code to find out.

[Andrew McRae]

"It shows up that this has serious consequences for form compilation" -- because the structure is lost, or something else?

[Jan Blechta]

@wence, it's not only about implementation details but rather about mathematics. Implementation details are consequences of it.

Suggested EnrichedElement is a nodal element built from nodal subelements. The enrichment (and the implementation in FIAT) is well-defined if and only if concatenated basis is linearly independent; see the test of ill-posed cases.

On the other hand suggested UnionElement just concatenates bases and does not produce meaningful nodes. FFC chain will not see any problem concatenating linearly dependent spaces. The problem will appear as singular system; try solving Poisson problem on P1+P1 (where + is current enriched element, i.e. UnionElement).

Summary: EnrichedElement is a nodal element of polynomial space being direct sum of given space (if well-defined). UnionElement is an element without nodes with polynomial space being union of given spaces. This is not an implementation detail; this is a clear mathematical definition.

[David Ham]

I don't see the advantage in gratuitously changing names. EnrichedElement has a current well-defined meaning in UFL. If you want a nodal enriched element, I suggest it would be better to create a new name for that.

[Lawrence Mitchell]

OK, but for the existing TensorProductElements that we build using EnrichedElement, these do have a unisolvent polynomial basis, it's just our implementation doesn't know how to make it right now. So mathematically EnrichedElement makes sense, but the implementation in this proposed scheme would require UnionElement.

IIUC, your FIAT pull request makes EnrichedElement do the following:

• orthogonalise the dual basis if the sub elements have a dual basis
• Throws an exception if the basis is not linearly independent.
• Falls back to basis concatenation if the sub elements don't have a dual basis implemented
• This may or may not a representation of a linearly independent basis (but you don't know).

Said change stops me from building a P1 + P0 element on triangles (say), because the basis is not linearly independent. Is your worry that you would like to continue to support that case?

[Jan Blechta]

I don't see the advantage in gratuitously changing names.

I'd like to have EnrichedElement and + produce nodal element in FEniCS as all other elements so far implemented.

EnrichedElement has a current well-defined meaning in UFL.

I would not agree here. UFL is not clear about the definition. Now it's hidden as an implementation detail that enriched is not nodal. I'd like to make this clear including the definition in EnrichedElement and UnionElement docstring.

[Jan Blechta]

but the implementation in this proposed scheme would require UnionElement

When you get the things right, you may switch back to EnrichedElement

This may or may not a representation of a linearly independent basis (but you don't know).

I do not understand this sentence.

You mean P1 + dP0. I haven't thought about this. How is it in fact implemented? Is one "extra" dof scrapped or not? Do you know a simple use, so that I can test this?

Is your worry that you would like to continue to support that case?

Primarily, the fallback to the old implementation of enriched was motivated by unit test of tensor product element.

[David Ham]

I definitely believe that + should not require a nodal basis. It just happens that all elements so far have nodal bases (or could have in the case of tensor product). UFL definitely should not outlaw or not support + on elements which are not nodal. For example trace elements are not nodal in the ordinary sense. Similarly, someone might choose to implement a modal finite element.

[Lawrence Mitchell]

This may or may not a representation of a linearly independent basis (but you don't know).

I do not understand this sentence.

Currently EnrichedElement doesn't make any guarantees on linear dependence (or not). As it happens, our uses with TensorProductElement always result in a linearly independent basis.

You mean P1 + dP0. I haven't thought about this. How is it in fact implemented? Is one "extra" dof scrapped or not? Do you know a simple use, so that I can test this?

I know some people seem to like P1 + dP0 for the pressure space in a P2-(P1 + dP0) discretisation for Stokes for mumble mumble local conservation mumble mumble reasons. But I can't find any of the references right now. It seems like a mad idea to me...

Primarily, the fallback to the old implementation of enriched was motivated by unit test of tensor product element.

As I understood it, making things work with this fallback didn't require any UFL changes, but maybe I'm missing something.

[Jan Blechta]

So that's another skeleton in the closet, whether DiscontinuousLagrangeTrace and HDivTrace are not nodal. They implement dual_basis().

a modal finite element

What is that?

[Lawrence Mitchell]

A reference for the stokes example I raised is, e.g.:

@Article{Boffi2012,
  author =       {Boffi, D.  and Cavallini, N.  and Gardini, F.  and
                  Gastaldi, L.},
  title =        {Local Mass Conservation of Stokes Finite Elements},
  journal =      {Journal of Scientific Computing},
  year =         2012,
  volume =       52,
  number =       2,
  pages =        {383--400},
  issn =         {1573-7691},
  doi =          {10.1007/s10915-011-9549-4},
  url =          {http://dx.doi.org/10.1007/s10915-011-9549-4}
}

In it they already note that although it seems to give better numerical properties the pressure mass matrix is rank deficient by 1 with this scheme: they solve using QR. The pressure laplacian is rank-deficient by the number of cells in the mesh...

[David Ham]

For modal finite element bases, see: http://onlinelibrary.wiley.com/doi/10.1002/nme.1620382204/full

[Jan Blechta]

I think that discussion shows that there is a need for the distinction in UFL. David objects strongly against renaming current EnrichedElement to UnionElement and changing the meaning of +. I'd be happy with keeping current EnrichedElement and '+', and introducing new NodalEnrichedElement. + could be deprecated later to avoid its ambiguity.

[David Ham]

I wonder if the ambiguity would be resolved more cleanly by turning this into a two operation process. So that EnrichedElement / + just concatenates basis functions, while Reorthogonalise(x + y) (or whichever name you like) would do the change of basis.

[Jan Blechta]

3d80006ff6c36ddd7a3912998ad02c9797b00713