Gauss-Legendre-Lobatto points used for the definition of the Lagrange elements
Issue #14
new
There is a general agreement that these interpolation points are better suited for high-order elements. It seems like this should be a direct modification of the Lagrange class (for someone who's familiar with FIAT.....).
Comments (9)
-
-
The issue is related to interpolation theory - Lagrange interpolants give different errors based on nodes used to construct them. A measure of this is the Lebesgue constant. You can see that equidistant nodes that Fiat uses are not the best choice. But of course, they are easy to implement. The problem with choosing different nodes is that in dims higher than one it is not known what the best choice is. Also these better nodes need to be computed. In 1d you have the luxury of Chebyshev points which are analytical and only slightly worse than the Gauss points. Higher dims are still researched, e.g Tim Warburton. Anyways I am not sure how big of a difference the nodes make for polynomial orders that are mostly used in fenics (which I assume is max 4). Finally @blechta GLL points are located on the boundaries so you'd enforce continuity in the usual way.
-
reporter
Thanks for following up on that, @blechta.
@mirok summed it up very well. I would just add that this can have an impact even with lower-order elements (2 or 3), e.g., when solving a wave equation propagating over a large number of wavelengths. Another advantage is that these points are also quadrature points. Using them for both interpolation and quadrature naturally leads to a diagonal mass matrix, hence avoiding lumping techniques.
-
Following @mirok's point about arbitrariness of nodes it seems to me that it would be quite useful allowing UFL and FFC chain to take user's subclass of
FIAT.FiniteElement. We have hit it already in #7 when different definition of dofs would be appreciated. Problem might be that FFC pretty much depends on a particular implementation. Opinions @martinal, @meg, @garth-wells? -
Note that you get good Lebesgue constants for higher degrees via Warburton's approach, but that for triangles & tets, the points are not also good quadrature points. The coordination of Lagrange nodes & quadrature points that the GLL trick gives is inherently one-dimensional or tensor-product. There is no lumping trick like that on the d-simplex for d > 1.
For degree 2 elements, the Warburton warp-blend nodes are the same as equispaced nodes. They differ starting at degree 3.
And also note that Warburton gets great numerical performance by doing things specialized for high-order DG that FEniCS and doesn't support (e.g. strong-form DG, pre-eliminating mass matrices). At orders for which FEniCS is going to be reasonably performant (low ones), there is a relatively small amount to be gained. That said, it would be a SMOP to dump other nodes into FIAT -- modepy by Andreas Warburton is open source, Python, and provides the warp-blend nodes.
-
By Andreas Warburton, Rob of course means Andreas Kloeckner. modeby is available here https://github.com/inducer/modepy.
-
Yes, Andreas Kloeckner. My bad.
-
Greetings all. The dham/gauss_lobatto branch contains a GLL implementation which I snarfed from Greg von Winckel. It's not currently in a state to merge but it does do most of the legwork on this issue. I had envisaged adding an additionall GLL element, but one could also just decide to replace the nodes in the Lagrange element. The branch also provides GLL quadrature, so one could apply the diagonal mass matrix trick (although FEniCS/Firedrake would not currently know that the matrix is diagonal so it would probably be more efficient to assemble v*dx and divide by that pointwise).
-
Pull request #14 provides GLL elements in 1D. The existing tensor product facilities therefore enable quad and hex elements based on them too.
- Log in to comment
Could you give more details? Are you thinking about continuous Lagrange elements? How does one ensure then interelement continuity?