Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #ifndef LIBMESH_FE_LAGRANGE_GLL_SHAPE_1D_H
20 : #define LIBMESH_FE_LAGRANGE_GLL_SHAPE_1D_H
21 :
22 : // Local includes
23 : #include "libmesh/libmesh_common.h"
24 : #include "libmesh/enum_order.h"
25 : #include "libmesh/int_range.h"
26 : #include "libmesh/quadrature_gauss_lobatto.h"
27 :
28 : // C++ includes
29 : #include <algorithm>
30 : #include <array>
31 : #include <iterator>
32 : #include <vector>
33 :
34 : // The one-dimensional nodal basis on the Gauss-Lobatto points. The tensor product element
35 : // types build their shape functions from it.
36 : //
37 : // The n points are -1, 1, and the roots of the derivative of the Legendre polynomial of
38 : // degree n-1, which is why they are also called the Gauss-Lobatto-Legendre (GLL) points.
39 : // The basis reads them from QGaussLobatto, so it and the rule with the same number of
40 : // points use identical values. A basis and a quadrature rule that share their points are
41 : // collocated: each shape function is one at its own point and zero at the others. The mass
42 : // matrix assembled with a collocated rule is therefore diagonal for any mapping and element
43 : // shape. That requires the points to match exactly, which is why the basis takes them from
44 : // the rule.
45 :
46 : namespace libMesh
47 : {
48 :
49 : /**
50 : * The one-dimensional nodal basis of degree n-1 on n Gauss-Lobatto points.
51 : */
52 : struct GLLBasis1D
53 : {
54 : /// The interpolation points in ascending order. The first is -1 and the last is 1.
55 : std::vector<Real> nodes;
56 :
57 : /// The barycentric weights. The i'th is the reciprocal of \f$\prod_{j \neq i}(\xi_i - \xi_j)\f$
58 : std::vector<Real> weights;
59 :
60 : /// \p deriv[k][i] is the derivative of the i'th shape function at the k'th point
61 : std::vector<std::vector<Real>> deriv;
62 :
63 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
64 : /// \p second[k][i] is the second derivative of the i'th shape function at the k'th point
65 : std::vector<std::vector<Real>> second;
66 : #endif
67 : };
68 :
69 :
70 : /**
71 : * The number of points in the largest Gauss-Lobatto rule that QGaussLobatto tabulates,
72 : * which it uses for orders FORTYSECOND and FORTYTHIRD. The highest available basis degree
73 : * is one less than this.
74 : */
75 : constexpr unsigned int fe_lagrange_gll_max_points_1D = 23;
76 :
77 :
78 : /**
79 : * \returns The one-dimensional nodal basis on \p n Gauss-Lobatto points.
80 : */
81 1210 : inline GLLBasis1D build_gll_basis_1D(const unsigned int n)
82 : {
83 44 : GLLBasis1D b;
84 :
85 : // QGaussLobatto chooses its number of points from the order: a rule of n points is exact
86 : // through degree 2n-3
87 1254 : QGaussLobatto rule(1, static_cast<Order>(2*n - 3));
88 1210 : rule.init(EDGE2);
89 44 : libmesh_assert_equal_to(rule.n_points(), n);
90 :
91 1210 : b.nodes.resize(n);
92 16335 : for (const auto i : make_range(n))
93 15675 : b.nodes[i] = rule.qp(i)(0);
94 :
95 1210 : b.weights.resize(n);
96 16335 : for (const auto i : make_range(n))
97 : {
98 550 : Real product = 1.;
99 252890 : for (const auto j : make_range(n))
100 237765 : if (j != i)
101 238832 : product *= b.nodes[i] - b.nodes[j];
102 15675 : b.weights[i] = 1. / product;
103 : }
104 :
105 : // Entry (k,i) is the derivative of shape function i at point k. Off the diagonal it is
106 : // (w_i/w_k)/(x_k - x_i). The shape functions sum to one, so each row sums to zero, and
107 : // the diagonal entry is minus the sum of the rest of its row. That also gives the
108 : // endpoint entries -p(p+1)/4 and p(p+1)/4, so they don't need a special case.
109 1254 : b.deriv.assign(n, std::vector<Real>(n, 0.));
110 16335 : for (const auto k : make_range(n))
111 : {
112 550 : Real row_sum = 0.;
113 252890 : for (const auto i : make_range(n))
114 237765 : if (i != k)
115 : {
116 263120 : b.deriv[k][i] = (b.weights[i] / b.weights[k]) / (b.nodes[k] - b.nodes[i]);
117 222640 : row_sum += b.deriv[k][i];
118 : }
119 16225 : b.deriv[k][k] = -row_sum;
120 : }
121 :
122 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
123 : // Off the diagonal, entry (k,i) of the second differentiation matrix is
124 : // 2 D_ki (D_kk - 1/(x_k - x_i)), with D the first differentiation matrix. Each row sums to
125 : // zero for the same reason as in D, which gives the diagonal.
126 1254 : b.second.assign(n, std::vector<Real>(n, 0.));
127 16335 : for (const auto k : make_range(n))
128 : {
129 550 : Real row_sum = 0.;
130 252890 : for (const auto i : make_range(n))
131 237765 : if (i != k)
132 : {
133 222640 : b.second[k][i] =
134 263120 : 2. * b.deriv[k][i] * (b.deriv[k][k] - 1. / (b.nodes[k] - b.nodes[i]));
135 222640 : row_sum += b.second[k][i];
136 : }
137 16225 : b.second[k][k] = -row_sum;
138 : }
139 : #endif
140 :
141 1254 : return b;
142 0 : }
143 :
144 : /**
145 : * \returns The one-dimensional basis of degree \p order.
146 : *
147 : * Every available basis is built in one static initialization on the first call. C++ makes
148 : * that initialization thread safe, so later calls read the bases without a lock. Each basis
149 : * of n points costs \f$O(n^2)\f$ operations, and building all of them takes well under a
150 : * millisecond.
151 : */
152 1006321128 : inline const GLLBasis1D & fe_lagrange_gll_1D_basis(const Order order)
153 : {
154 55 : static const std::array<GLLBasis1D, fe_lagrange_gll_max_points_1D + 1> bases = []()
155 : {
156 53 : std::array<GLLBasis1D, fe_lagrange_gll_max_points_1D + 1> b;
157 1265 : for (const auto n : make_range(2u, fe_lagrange_gll_max_points_1D + 1))
158 1210 : b[n] = build_gll_basis_1D(n);
159 55 : return b;
160 1006321128 : }();
161 :
162 1006321128 : const unsigned int n = static_cast<unsigned int>(order) + 1;
163 :
164 1006321128 : libmesh_error_msg_if(n < 2 || n > fe_lagrange_gll_max_points_1D,
165 : "The Gauss-Lobatto nodal basis supports orders 1 through "
166 : << fe_lagrange_gll_max_points_1D - 1 << ", but order " << order
167 : << " was requested.");
168 :
169 1090181222 : return bases[n];
170 : }
171 :
172 :
173 : /**
174 : * \returns The index of the point of \p b equal to \p xi, or \p invalid_uint if there is
175 : * none.
176 : *
177 : * The points are in ascending order, so a binary search finds the candidate. The comparison
178 : * is exact. The points of a collocated quadrature rule still match, because the basis and
179 : * the rule read the same values from QGaussLobatto.
180 : */
181 1006321128 : inline unsigned int fe_lagrange_gll_1D_point(const GLLBasis1D & b, const Real xi)
182 : {
183 1006321128 : const auto it = std::lower_bound(b.nodes.begin(), b.nodes.end(), xi);
184 :
185 1006321128 : if (it != b.nodes.end() && *it == xi)
186 115896480 : return cast_int<unsigned int>(std::distance(b.nodes.begin(), it));
187 :
188 74202054 : return invalid_uint;
189 : }
190 :
191 : /**
192 : * \returns The value of the \p i'th one-dimensional Gauss-Lobatto nodal shape function of
193 : * degree \p order at \p xi.
194 : *
195 : * At a basis point the value is exactly one or zero, which a collocated quadrature rule
196 : * relies on. Elsewhere it is evaluated with the barycentric formula
197 : * \f$\ell_i(\xi) = \frac{w_i/(\xi - \xi_i)}{\sum_j w_j/(\xi - \xi_j)}\f$.
198 : */
199 619912104 : inline Real fe_lagrange_gll_1D_shape(const Order order,
200 : const unsigned int i,
201 : const Real xi)
202 : {
203 619912104 : const GLLBasis1D & b = fe_lagrange_gll_1D_basis(order);
204 :
205 51659342 : libmesh_assert_less (i, b.nodes.size());
206 :
207 619912104 : const unsigned int k = fe_lagrange_gll_1D_point(b, xi);
208 619912104 : if (k != invalid_uint)
209 146879272 : return (i == k) ? 1. : 0.;
210 :
211 44505766 : Real numerator = 0., denominator = 0.;
212 :
213 2875495680 : for (const auto j : index_range(b.nodes))
214 : {
215 2536545362 : const Real term = b.weights[j] / (xi - b.nodes[j]);
216 :
217 2341426488 : denominator += term;
218 :
219 2341426488 : if (j == i)
220 44505766 : numerator = term;
221 : }
222 :
223 534069192 : return numerator / denominator;
224 : }
225 :
226 : /**
227 : * \returns The barycentric interpolant at \p xi of the values \p f[k] at the points of
228 : * \p b, for \p xi away from those points.
229 : *
230 : * This is \f$\sum_k \ell_k(\xi) f_k\f$ in the second barycentric form, which stays
231 : * accurate as \p xi approaches a point.
232 : */
233 : template <typename Column>
234 356355456 : inline Real fe_lagrange_gll_1D_interpolate(const GLLBasis1D & b,
235 : const Real xi,
236 : const Column & f)
237 : {
238 29696288 : Real numerator = 0., denominator = 0.;
239 :
240 1918430880 : for (const auto k : index_range(b.nodes))
241 : {
242 1692248376 : const Real term = b.weights[k] / (xi - b.nodes[k]);
243 :
244 1692248376 : numerator += term * f(k);
245 1562075424 : denominator += term;
246 : }
247 :
248 356355456 : return numerator / denominator;
249 : }
250 :
251 : /**
252 : * \returns The derivative of the \p i'th one-dimensional Gauss-Lobatto nodal shape function
253 : * of degree \p order at \p xi.
254 : *
255 : * At a basis point it is read from the differentiation matrix. The derivative is a
256 : * polynomial of degree p-1, which the basis represents exactly, so elsewhere it is the
257 : * interpolant of column \p i of that matrix. Differentiating the barycentric formula
258 : * directly would instead cancel catastrophically as \p xi approaches a point.
259 : */
260 265169388 : inline Real fe_lagrange_gll_1D_shape_deriv(const Order order,
261 : const unsigned int i,
262 : const unsigned int libmesh_dbg_var(j),
263 : const Real xi)
264 : {
265 : // only d()/dxi in 1D!
266 24106308 : libmesh_assert_equal_to (j, 0);
267 :
268 265169388 : const GLLBasis1D & b = fe_lagrange_gll_1D_basis(order);
269 :
270 24106308 : libmesh_assert_less (i, b.nodes.size());
271 :
272 265169388 : const unsigned int k = fe_lagrange_gll_1D_point(b, xi);
273 265169388 : if (k != invalid_uint)
274 20619126 : return b.deriv[k][i];
275 :
276 266782104 : return fe_lagrange_gll_1D_interpolate(b, xi, [&b, i](const std::size_t m)
277 1386689454 : { return b.deriv[m][i]; });
278 : }
279 :
280 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
281 : /**
282 : * \returns The second derivative of the \p i'th one-dimensional Gauss-Lobatto nodal shape
283 : * function of degree \p order at \p xi.
284 : *
285 : * At a basis point it is read from the second differentiation matrix. Elsewhere it is the
286 : * interpolant of column \p i of that matrix, for the same reason as the first derivative.
287 : */
288 89038884 : inline Real fe_lagrange_gll_1D_shape_second_deriv(const Order order,
289 : const unsigned int i,
290 : const unsigned int libmesh_dbg_var(j),
291 : const Real xi)
292 : {
293 : // only d^2()/dxi^2 in 1D!
294 8094444 : libmesh_assert_equal_to (j, 0);
295 :
296 89038884 : const GLLBasis1D & b = fe_lagrange_gll_1D_basis(order);
297 :
298 8094444 : libmesh_assert_less (i, b.nodes.size());
299 :
300 89038884 : const unsigned int k = fe_lagrange_gll_1D_point(b, xi);
301 89038884 : if (k != invalid_uint)
302 6929978 : return b.second[k][i];
303 :
304 89573352 : return fe_lagrange_gll_1D_interpolate(b, xi, [&b, i](const std::size_t m)
305 465428162 : { return b.second[m][i]; });
306 : }
307 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
308 :
309 : } // namespace libMesh
310 :
311 : #endif // LIBMESH_FE_LAGRANGE_GLL_SHAPE_1D_H
|
