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_H
20 : #define LIBMESH_FE_LAGRANGE_GLL_SHAPE_H
21 :
22 : // Local includes
23 : #include "libmesh/enum_fe_family.h"
24 : #include "libmesh/enum_order.h"
25 : #include "libmesh/fe_lagrange_gll_shape_1D.h"
26 : #include "libmesh/int_range.h"
27 :
28 : // Tensor product evaluation shared by the Gauss-Lobatto nodal families. The families differ
29 : // only in how a shape function index maps to a tensor index, so each family specializes
30 : // fe_lagrange_gll_tensor_index and shares the rest.
31 :
32 : namespace libMesh
33 : {
34 :
35 : /**
36 : * The pair of coordinates differentiated by each second derivative index \p j of
37 : * FE::shape_second_deriv().
38 : *
39 : * In dimension d only the first \f$d(d+1)/2\f$ entries are used: one on an edge, three on
40 : * a face, and all six on a cell.
41 : */
42 : constexpr unsigned int fe_lagrange_gll_second_deriv_pairs[6][2] =
43 : {{0, 0}, // d^2 phi / dxi^2
44 : {0, 1}, // d^2 phi / dxi deta
45 : {1, 1}, // d^2 phi / deta^2
46 : {0, 2}, // d^2 phi / dxi dzeta
47 : {1, 2}, // d^2 phi / deta dzeta
48 : {2, 2}}; // d^2 phi / dzeta^2
49 :
50 :
51 : /**
52 : * \returns The number of distinct second derivatives in dimension \p Dim.
53 : */
54 : template <unsigned int Dim>
55 16100376 : constexpr unsigned int fe_lagrange_gll_n_second_derivs()
56 : {
57 16100376 : return Dim * (Dim + 1) / 2;
58 : }
59 :
60 : /**
61 : * Writes the tensor index of shape function \p i of order \p order into the first \p dim
62 : * entries of \p indices. Entry d is the index of the one-dimensional shape function in
63 : * coordinate d.
64 : *
65 : * Each Gauss-Lobatto nodal family specializes this function. Using a family without a
66 : * specialization is a link error.
67 : */
68 : template <FEFamily T>
69 : void fe_lagrange_gll_tensor_index(const Elem * elem,
70 : const Order order,
71 : const unsigned int dim,
72 : const unsigned int i,
73 : unsigned int * indices);
74 :
75 :
76 : /**
77 : * All L2_LAGRANGE_GLL degrees of freedom belong to the element. They are numbered
78 : * lexicographically over the tensor grid, with the first coordinate varying fastest. QBase
79 : * orders the points of a tensor product rule the same way, so shape function i is the one
80 : * that is one at quadrature point i. The numbering doesn't depend on the element, so the
81 : * element argument is unused.
82 : */
83 : template <>
84 28167913 : inline void fe_lagrange_gll_tensor_index<L2_LAGRANGE_GLL>(const Elem *,
85 : const Order order,
86 : const unsigned int dim,
87 : const unsigned int i,
88 : unsigned int * indices)
89 : {
90 309847043 : const unsigned int n = static_cast<unsigned int>(order) + 1;
91 :
92 28167913 : unsigned int stride = 1;
93 :
94 1344336084 : for (const auto d : make_range(dim))
95 : {
96 1006321128 : indices[d] = (i / stride) % n;
97 1006321128 : stride *= n;
98 : }
99 :
100 28167913 : libmesh_assert_less (i, stride);
101 28167913 : }
102 :
103 : /**
104 : * \returns Shape function \p i of family \p T in dimension \p Dim, differentiated once
105 : * with respect to each coordinate listed in the first \p n_deriv_coords entries of
106 : * \p deriv_coords.
107 : *
108 : * Each factor of the tensor product is a one-dimensional shape function, differentiated as
109 : * many times as its coordinate is listed. FE::shape() passes zero coordinates,
110 : * FE::shape_deriv() passes one, and FE::shape_second_deriv() passes two.
111 : */
112 : template <FEFamily T, unsigned int Dim>
113 338014956 : inline Real fe_lagrange_gll_tensor_product(const Elem * elem,
114 : const Order order,
115 : const unsigned int i,
116 : const Point & p,
117 : const unsigned int * const deriv_coords,
118 : const unsigned int n_deriv_coords)
119 : {
120 : unsigned int indices[3];
121 :
122 28167913 : fe_lagrange_gll_tensor_index<T>(elem, order, Dim, i, indices);
123 :
124 28167913 : Real value = 1.;
125 :
126 1344336084 : for (const auto d : make_range(Dim))
127 : {
128 83860094 : unsigned int deriv_order = 0;
129 2448521160 : for (const auto k : make_range(n_deriv_coords))
130 1442200032 : if (deriv_coords[k] == d)
131 483542352 : ++deriv_order;
132 :
133 1006321128 : switch (deriv_order)
134 : {
135 619912104 : case 0:
136 619912104 : value *= fe_lagrange_gll_1D_shape(order, indices[d], p(d));
137 619912104 : break;
138 :
139 289275696 : case 1:
140 289275696 : value *= fe_lagrange_gll_1D_shape_deriv(order, indices[d], 0, p(d));
141 289275696 : break;
142 :
143 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
144 97133328 : case 2:
145 97133328 : value *= fe_lagrange_gll_1D_shape_second_deriv(order, indices[d], 0, p(d));
146 97133328 : break;
147 : #endif
148 :
149 0 : default:
150 0 : libmesh_error_msg("Invalid derivative order " << deriv_order << ".");
151 : }
152 : }
153 :
154 338014956 : return value;
155 : }
156 :
157 : } // namespace libMesh
158 :
159 : #endif // LIBMESH_FE_LAGRANGE_GLL_SHAPE_H
|
