Repository navigation
Add the L2_LAGRANGE_GLL family - #4583
Conversation
c2a37aa to
658d00a
Compare
The EDGE4 case subtracted a Point from nodes[3] and discarded the result, so the node stayed at the origin. It now sits at xi = 1/3, matching Edge4::master_point(). Co-Authored-By: Claude <noreply@anthropic.com>
The basis of degree p interpolates at the p+1 Gauss-Lobatto points: -1, 1, and the roots of the derivative of the Legendre polynomial of degree p, also called the Gauss-Lobatto-Legendre points. It reads them from QGaussLobatto, so it shares its points exactly with the Gauss-Lobatto rule of the same size. At a basis point, a shape function is exactly one or zero, and its derivatives are read from precomputed differentiation matrices. The point lookup is a binary search on the ascending points. A collocated quadrature rule therefore gets its shape function values and derivatives without any polynomial evaluation. At other points, the value comes from the barycentric formula, and each derivative is the barycentric interpolant of a column of its differentiation matrix. That is exact, because the basis represents a derivative of degree p-1 or p-2 exactly, and it stays accurate as the point approaches a basis point, where differentiating the barycentric formula directly would cancel catastrophically. Off-diagonal entries of the differentiation matrix are ratios of barycentric weights over point differences, and those of the second differentiation matrix follow from them by a closed form. The shape functions sum to one, so each row of either matrix sums to zero, and each diagonal entry is minus the sum of the rest of its row. That also produces the endpoint entries -p(p+1)/4 and p(p+1)/4 of the first. Building a basis of n points costs O(n^2) operations. The second differentiation matrix is built only when LIBMESH_ENABLE_SECOND_DERIVATIVES is set. All bases are built in one function-local static on first use, which C++ initializes thread safely, so later lookups need no lock. Co-Authored-By: Claude <noreply@anthropic.com>
L2_LAGRANGE_GLL is a discontinuous nodal basis on the Gauss-Lobatto points. Each shape function is one at its own point and zero at the others, so the mass matrix assembled with the Gauss-Lobatto rule on those points is diagonal for any mapping and element shape. The basis is a tensor product, so it is defined only on tensor product element types. All of its degrees of freedom belong to the element. Above order two the Gauss-Lobatto points are not element nodes, so the order is independent of the node count; order eight works on a QUAD4, for example. Supported orders run from one to one less than the number of points in the largest tabulated Gauss-Lobatto rule, and higher orders raise an error. FEType::default_quadrature_rule() returns the collocated rule. A Gauss-Lobatto rule of n points is exact through degree 2n-3, so a basis of degree p, with p+1 points, needs order 2p-1. This also holds under p refinement: QBase raises a rule's order by twice the element's p_level, and 2p-1 rises by two per degree. The collocated rule is exact through degree 2p-1 in each coordinate. The mass integrand has degree 2p in each coordinate, so it is underintegrated, which is what makes the mass matrix diagonal. In one dimension the stiffness integrand has degree 2p-2 and is integrated exactly on an affine element; in two and three dimensions its undifferentiated factors still have degree 2p, so it is underintegrated in the same way. Neither costs the optimal convergence rates on affine meshes, O(h^p) in H1 and O(h^(p+1)) in L2, which is the standard result for Gauss-Lobatto collocation in spectral element methods. Co-Authored-By: Claude <noreply@anthropic.com>
The existing FETest harness covers the shape functions, their first and second derivatives, and the degree of freedom counts, at orders one through four on every tensor product element type, serendipity types included. All degrees of freedom belong to the element, so every order works on every one of those types. The new FEGLLCollocationTest checks the collocation property. At the points of the family's default quadrature rule, the matrix of shape function values is the identity, and the assembled mass matrix is diagonal with entries equal to JxW. Because the basis reads its points from QGaussLobatto, these comparisons are exact and use no tolerance. The test mesh is sheared so that each element's mapping has off-diagonal Jacobian terms. testPartitionOfUnity requires every family to be listed explicitly, so this adds L2_LAGRANGE_GLL to the families whose shape functions sum to one. Co-Authored-By: Claude <noreply@anthropic.com>
The tensor product evaluation and the table of second derivative coordinate pairs were duplicated in each of the three dimension files, and the continuous LAGRANGE_GLL family in libMesh#4552 would have duplicated them again. They now live once in fe_lagrange_gll_shape.h, as fe_lagrange_gll_tensor_product(). The dimension files keep only the per-dimension FE specializations, as fe_lagrange_shape_2D.C does for LAGRANGE and L2_LAGRANGE. The families differ only in how a shape function index maps to a tensor index, so that mapping is a function template, fe_lagrange_gll_tensor_index, that each family specializes. Using a family without a specialization is a link error. A single second derivative pair table replaces the per-dimension tables. In dimension d, FE::shape_second_deriv() uses the first d(d+1)/2 indices, so an edge uses one entry, a face three, and a cell all six. Co-Authored-By: Claude <noreply@anthropic.com>
At orders one and two on the tensor product element types, the Gauss-Lobatto
points {-1,1} and {-1,0,1} are the element nodes, so the Gauss-Lobatto nodal
families and LAGRANGE have the same shape functions. The implementations are independent:
LAGRANGE evaluates a closed form for each order, and the Gauss-Lobatto basis
evaluates the barycentric formula on its point table.
The new test pairs each Gauss-Lobatto shape function with the LAGRANGE shape
function of the node at the same reference point and checks that their values agree on a
grid of sample points. It catches errors in the point table, the barycentric evaluation,
and the tensor index map.
QUAD8 and HEX20 are excluded. On those serendipity types LAGRANGE has 8 and 20 shape
functions, against 9 and 27 for the tensor product basis of the same order.
Co-Authored-By: Claude <noreply@anthropic.com>
roystgnr
left a comment
There was a problem hiding this comment.
Let's wait until Monday afternoon to merge, on the off chance @jwpeterson wants to review and/or has any better ideas for the enum name.
Since this is ready to merge there's no need to split out the EDGE4 refspace fix, but we still might want to sic Claude on figuring out why @rochi00 's recent test coverage still missed it. I recall thinking it looked pretty thorough!
| using namespace libMesh; | ||
|
|
||
| template <ElemType elem_type, Order order> | ||
| class FEGLLCollocationTest : public CppUnit::TestCase |
There was a problem hiding this comment.
I really like this - would it be reasonable to generalize it a little and instantiate it for LAGRANGE/L2_LAGRANGE and QNODAL as well?
| {"LAGRANGE_VEC" , LAGRANGE_VEC}, | ||
| {"L2_LAGRANGE" , L2_LAGRANGE}, | ||
| {"L2_LAGRANGE_VEC" , L2_LAGRANGE_VEC}, | ||
| {"L2_LAGRANGE_GLL" , L2_LAGRANGE_GLL}, |
There was a problem hiding this comment.
Any chance this is going to confuse people? Yes, this FE has the Lagrange property in the sense of "its canonical basis is given by 'delta function' functionals at a set of points", and in the sense of "using those point values as coefficients gives you the Lagrange interpolating polynomial"... but I worry that too many MOOSE users might think "Lagrange" specifically means "interpolates at the nodes".
I can't think of an obviously better name, though. L2_CARDINAL_GLL works, since "φᵢ(xⱼ)=δᵢⱼ" also gets called the "cardinal property", but that's much much less common a name AFAIK. The only thing going for it is that people will correctly think "I should look up what that means" rather than incorrectly thinking "I already understand what that means".
There was a problem hiding this comment.
"I already understand what that means".
I don't think there's much danger of MOOSE users thinking this. I guess I can only hypothesize but I feel like a run-of-the-mill MOOSE user who is used to interacting with LAGRANGE, if they were to see LAGRANGE_GLL would not think they understand what it means. When I first saw Gauss-Lobatto in MFEM/libMesh contexts I did not understand what it meant, e.g. I had to learn that Gauss implicitly means Gauss-Legendre with the points as roots of the
Now you're right in that I did think I knew what Lagrange meant ... until I worked on this collocated GLL space and learned it also fit the definition of Lagrange. I mean it seems obviously in hind-sight that a basis should be defined by mathematical properties and not an implementation detail that the interpolation points happen to be at things called Nodes.
Anyway it's my opinion that the name is a good thing in that it will challenge users preconception of what it means to be Lagrange; I think a challenged preconception is much more likely than the an belief in understanding.
| const unsigned int i, | ||
| const Point & p) | ||
| { | ||
| return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1> |
There was a problem hiding this comment.
Would it be safer to just call e.g. fe_lagrange_gll_1D_shape in the 1D cases? I'm 99% sure the compiler will optimize away all the unnecessary loops but it wouldn't hurt to get to 100%.
| nodes[0] = Point (-1.,0.,0.); | ||
| nodes[1] = Point (1.,0.,0.); | ||
| nodes[2] = Point (-1./3.,0.,0.); | ||
| nodes[3] - Point (1./3.,0.,0.); |
There was a problem hiding this comment.
What. The. Hell.
This would be worthy of a PR on its own. And didn't Rochi just finally add a bunch of refspace-nodes test coverage for us after discovering issues there? How did we miss this one?
Adds
L2_LAGRANGE_GLL, a discontinuous nodal basis on the Gauss-Lobatto points of the tensor-product element types. Its shape functions are collocated with the Gauss-Lobatto rule of the same order, so that rule assembles an exactly diagonal mass matrix for any mapping. All DOFs belong to the element, so the order is independent of the node count (e.g. order 8 on a QUAD4).FEType::default_quadrature_rule()returns the collocated rule. Under that rule stiffness stays exact on affine elements, while mass is underintegrated by one degree; that is what makes it diagonal.The 1D basis takes its points from a one-dimensional
QGaussLobattorule, so the basis and the collocated rule carry the same values. Tests cover the shape functions through second derivatives, the collocation property, and agreement withLAGRANGEat orders 1 and 2, where the two bases coincide. The tensor-product evaluation is shared infe_lagrange_gll_shape.hso the continuous family in #4552 can reuse it.🤖 Generated with Claude Code