Better nodes for RT/Nedelec

[wence-]

Background

This is related to #15 and also https://bitbucket.org/fenics-project/fiat/issues/25/nedelec-dofs-give-suboptimal-interpolation

Although the nodes that FIAT implements for these elements give optimal order of convergence for projection, they are not sufficient to provide optimal interpolation order. For example, for RTq, the textbook nodes are:

1. \int_F v \cdot n p ds where p is polynomials of degree q-1;
2. \int_T v \cdot p dx where p is d-vector-valued polynomials of degree q-2.

With these nodes, the so defined interpolation operator I^q_T has the approximation property (e.g. Brezzi and Fortin, Chapter III.3 (1991)):

||u - I^q_T u||_H(div, T) <= C h^q_T |u|_H^{q+1}_T

||u - I^q_T u||L^2(T) <= C h^q_T |u|_H^q_T

tl;dr: the interpolation error should converge at order q in both the L2 and Hdiv norms.

Instead of integral moments, FIAT uses point normal evaluations on the faces and point component evaluations on the cell (at equispaced lattice points). We therefore lose one order of approximation in interpolation in the Hdiv norm (because effectively we're underintegrating the moments).

On triangles, the face (edge) dofs can be component evaluations if they are chosen to fall at the Gauss points (and the polynomial set is the Gauss-Legendre polynomials) since then the quadrature rule collocates with the nodes of polynomial. No such construction exists for face dofs on tetrahedra.

Proposed fix

To fix this, one should replace the underintegrated point evaluation nodes by textbook integral moments. I've done this in firedrakeproject@c2e2eeb and obtain correct interpolation order for RTq.

Policy implications

Changing FIAT's dual basis modifies the primal basis, which means any checkpoints a user has now don't work. We therefore can't make this change without providing a migration path, and also deprecation warnings.

Proposal

UFL finite elements already have variant tags so that one can (for example) select between equispaced and Gauss-Lobatto points on interval elements. We should introduce variants (say integral [new] and pointwise [old]) defaulting to pointwise for the elements we fix. The form compiler's UFL->fiat element translation then needs to handle this. FIAT should spew deprecation warnings for the pointwise elements and suggest the fix for migrating data (project from old definition to new definition). At some point later (next major fenics release after the deprecation warning is introduced?) we switch the default. In the following release we remove the old nodes.

[dham]

Agree. I suggest that on the Firedrake side we use our usual more aggressive strategy of giving month's warning and then changing the default on master.

[jhale]

Fully supportive of the change.

This was discussed in 2015 with @rckirby in the context of projection (and hopefully interpolation if it was done correctly!) operators for the Reissner-Mindlin plate problem see https://www.sciencedirect.com/science/article/abs/pii/0045794989903659

I hope Rob won't mind if I share the emails here as I think they are useful for the discussion:

Hi Rob,

I'm trying to understand something in FIAT and connect it with the
definition in the FEniCS book. The reason I need to understand is
related to implementing MITC plate elements in FEniCS, which involves
construction of a local projection/tying of a vector-valued CG space
to a NED space. I was hoping you could shed some light on it as you
wrote the code and are an expert in such matters.

In the FEniCS book, the edge DOFs are defined as:

\int_{e} v \cdot t p for p \in P_{q-1)(e) for each edge e

So, if we focus in on NED_2^1 I would expect to be able to
equivalently replace this integral with point evaluation at two Gauss
quadrature points on the edge.

However, my reading of FIAT source code:

https://bitbucket.org/fenics-project/fiat/src/81e0314b92c63b44d0ad06918f6d1eafa9528b23/FIAT/nedelec.py?at=master#cl-175

is that the point evaluations are evaluated using points generated
using make_points, which are uniformly distributed along the edge.

Whereas for the internal moments, the evaluations are indeed done
using gauss quadrature points.

https://bitbucket.org/fenics-project/fiat/src/81e0314b92c63b44d0ad06918f6d1eafa9528b23/FIAT/nedelec.py?at=master#cl-184

Why is it done this way? What are the differences between the two
approaches? Perhaps we've missed something obvious?

Thanks for your help,
Jack S. Hale

Jack,
Thanks for the question.
It's been a while since I've dug into this.
1.) Any set of (the right number of unisolvent) points on the edge is
fine, as well the moments. You would gain some stability advantages at
high order and maybe some cute discrete orthogonality if you use Gauss
points on the edges. There isn't a comparable "magic" set of points for
2d facets of 3d elements, so I didn't put too much work into optimizing
it in the 2d case. Point values vs moments matter some in how the
elements transform under pullback (you need to account for face Jacobians).
2.) I'm using gaussian integration for the internal moments because
those are integrals, and equispaced points make a bad integration rule.

Does this make sense?
Rob

Actually, the 'cute discrete orthogonality' was exactly what we were
looking for I think. We are trying to construct an element-level
projection operator relating the DOFs of [CG_2]^2+[B_3]^2 to NED_2^1.
As a first check, we were trying to construct a projection of NED_2^1
onto itself. So we were expecting to get an identity matrix. But we
did not. This is because we were trying to construct the tying on the
assumption that the point evaluations were at the Gauss quadrature
points. This is what is typically done in the 'MITC literature'.

We have been successful using this approach to construct an operator
from [CG_2]^2 to NED_1^1 and making a Duran-Liberman type finite
element for simulating Reissner-Mindlin plates. This is because the
Gauss Quadrature rule and the point evaluations coincide for this
simple case.

Jack

[rckirby]

I am simply mortified that Jack would reveal private correspondence! ;) But, I am even more mortified (mea maxima culpa!), that perhaps I didn't appreciate until recently that there are actually two issues at stake in how the boundary degrees of freedom until this morning: 1.) Ensuring enough continuity between cells that piecewise polynomials are actually in H(div) or H(curl). That plus having the right polynomials in the space lets you get optimal variational problems. 2.) Determining the accuracy of the Ciarlet-type nodal interpolation operator defined by those DOFs. Note that for full accuracy in variational problems, you don't need to have DOFs that provide optimality, you just need such a set of DOFs to exist. IIRC, in neither my FIAT papers nor the FFC paper with Marie on div/curl-conforming elements did we actually investigate the accuracy of the nodal interpolant, just the result of solving finite element variational problems. Now, I claim that the "dumb" nodes currently in FIAT for Raviart-Thomas/Nedelec are actually fine for solving variational problems. However, it seems very likely that if I have a given smooth vector field u, then knowing (u.n) at given points on an edge \gamma may not specify \int_\gamma (u.n) \mu ds on an edge for polynomials \mu of high enough degree. In fact, it will be right for RT elements iff we choose the Gauss points. The exact error estimate one can expect for the interpolant in the absence of full quadrature accuracy of the points is probably somewhat subtle. So, I also support changing the boundary degrees of freedom to use moments against polynomials of a certain degree (you can hack the 2d case by switching to Gauss points, but then you are sunk in 3d). The infrastructure exists in FIAT to do this pretty simply (in fact, others have branches doing most of it already), although there will be a slight slow-down in instantiating elements. FYI, I have a languishing FIAT branch that should speed this up, but last time we tried to merge some time ago it broke a few things. This should be revisited when I get a chance, but doesn't affect correctness of switching to moments, just performance. Best regards, Rob Kirby

…

On Thu, Apr 2, 2020 at 6:43 AM Jack S. Hale ***@***.***> wrote: Fully supportive of the change. This was discussed in 2015 with @rckirby <https://github.com/rckirby> in the context of projection (and hopefully interpolation if it was done correctly!) operators for the Reissner-Mindlin plate problem c.f. https://www.sciencedirect.com/science/article/abs/pii/0045794989903659 I hope Rob won't mind if I share the emails here as I think they are useful for the discussion: Hi Rob, I'm trying to understand something in FIAT and connect it with the definition in the FEniCS book. The reason I need to understand is related to implementing MITC plate elements in FEniCS, which involves construction of a local projection/tying of a vector-valued CG space to a NED space. I was hoping you could shed some light on it as you wrote the code and are an expert in such matters. In the FEniCS book, the edge DOFs are defined as: \int_{e} v \cdot t p for p \in P_{q-1)(e) for each edge e So, if we focus in on NED_2^1 I would expect to be able to equivalently replace this integral with point evaluation at two Gauss quadrature points on the edge. However, my reading of FIAT source code: https://bitbucket.org/fenics-project/fiat/src/81e0314b92c63b44d0ad06918f6d1eafa9528b23/FIAT/nedelec.py?at=master#cl-175 is that the point evaluations are evaluated using points generated using make_points, which are uniformly distributed along the edge. Whereas for the internal moments, the evaluations are indeed done using gauss quadrature points. https://bitbucket.org/fenics-project/fiat/src/81e0314b92c63b44d0ad06918f6d1eafa9528b23/FIAT/nedelec.py?at=master#cl-184 Why is it done this way? What are the differences between the two approaches? Perhaps we've missed something obvious? Thanks for your help, Jack S. Hale Jack, Thanks for the question. It's been a while since I've dug into this. 1.) Any set of (the right number of unisolvent) points on the edge is fine, as well the moments. You would gain some stability advantages at high order and maybe some cute discrete orthogonality if you use Gauss points on the edges. There isn't a comparable "magic" set of points for 2d facets of 3d elements, so I didn't put too much work into optimizing it in the 2d case. Point values vs moments matter some in how the elements transform under pullback (you need to account for face Jacobians). 2.) I'm using gaussian integration for the internal moments because those are integrals, and equispaced points make a bad integration rule. Does this make sense? Rob Actually, the 'cute discrete orthogonality' was exactly what we were looking for I think. We are trying to construct an element-level projection operator relating the DOFs of [CG_2]^2+[B_3]^2 to NED_2^1. As a first check, we were trying to construct a projection of NED_2^1 onto itself. So we were expecting to get an identity matrix. But we did not. This is because we were trying to construct the tying on the assumption that the point evaluations were at the Gauss quadrature points. This is what is typically done in the 'MITC literature'. We have been successful using this approach to construct an operator from [CG_2]^2 to NED_1^1 and making a Duran-Liberman type finite element for simulating Reissner-Mindlin plates. This is because the Gauss Quadrature rule and the point evaluations coincide for this simple case. Jack — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#40 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA4XUL4ZB6UNC7RLWUX2TZDRKR25HANCNFSM4L2ID4YA> .

[jhale]

Don't worry Rob I'll keep the other things under lock and key ;o)

Until @wence- mentioned it today I was also not aware of the effect of this choice on the convergence of the global interpolant. Interesting to note that we had no convergence issues when implementing MITC7 which uses the local interpolant in the context of assembling a global variational problem.

[pefarrell]

I agree that we should use the better nodes. I support the transition plan proposed by Lawrence, and the more aggressive strategy for firedrake proposed by David. What are the next steps?

[wence-]

PR opened (#52).