Definition of degrees of freedom for Nedelec element of first kind

[Jack Hale]

@meg @robert_kirby @cmaurini @lzlarryli @spike05

This report is based around a private conversation I had with Rob Kirby. I thought it might be worth discussing it here as well, and whether changes should be made to FIAT.

I will just talk about the family NED^1 on triangles but it may well apply elsewhere in FIAT.

At the moment, the edge degrees of freedom for NED_^ are located at equally spaced points along the edge of the triangle. So for NED_2^1 we have points at (1/3, 2/3) that are passed to PointEdgeTangentEvaluation to make NedelecDual2D.

As Rob pointed out quite rightly, these are a unisolvent set of points, and it works fine.

However, it is also possible to define the edge DOFs of NED_2^1 as an evaluation at the quadrature points of a Gauss rule of order 2. In my opinion, this definition is more consistent with the integral definition of NED^1 family as given in the FEniCS book.

The practical reason for doing this is that we then have a nice Galerkin-orthogonality property which is very useful for defining the local projection of the degrees of freedom of a vector-valued CG function onto a NED^1 space. By a fluke we were successful doing this for [CG_2]^2 to NED^1_1 space, because the FIAT definition and the integral definition coincide. These projections are critical to the success of our fenics-shells project.

Rob also mentioned that conditioning number for higher-order NED^1 elements may improve as well.

Proposal:

1. Make a class in functional.py IntegralEdgeTangentMomentEvaluation
2. Adjust NedelecDual2D definition appropriately.

We can probably do the coding on this but would need a bit of guidance.

[Jack Hale]

Flipped family and order consistent with book.

[Marie Elisabeth Rognes]

I don't have a strong opinion on the choice of degrees of freedom, so if there is a better choice, please go ahead. The only thing I would keep in mind is to check the resulting performance for the generated code in the relevant places. For higher orders, the differences between a single point evaluation and a possible large number in a high order quadrature rule may be significant.

[Jan Blechta]

Or there could be a user interface for this (while not modifying current definition)

fiat_element = FIAT.FiniteElement(...)
ufl_element = ufl.FiniteElement(fiat_element)
V = FunctionSpaceBase(mesh, ufl_element)

This would need a work in FFC and UFL. But maybe it is not correct abstraction as UFL elements do not know anything about FIAT elements now.