libMesh/libmesh: coverage diff

Base 9e5b0b Head #4583 689827
Total Total +/- New
Rate 66.51% 66.51% +0.00% 100.00%
Hits 80359 80584 +225 337
Misses 40472 40583 +111 0
Filename Stmts Miss Cover
include/fe/fe_lagrange_gll_shape.h +31 +2 +93.55%
include/fe/fe_lagrange_gll_shape_1D.h +87 +1 +98.85%
include/fe/fe_lagrange_shape_1D.h 0 -14 +13.21%
src/fe/fe_abstract.C +8 -4 +1.03%
src/fe/fe_base.C +8 +2 +0.11%
src/fe/fe_boundary.C +1 +1 -0.18%
src/fe/fe_interface.C +6 +2 +0.12%
src/fe/fe_l2_lagrange_gll.C +77 +42 +45.45%
src/fe/fe_lagrange_gll_shape_0D.C +22 +22 +100.00%
src/fe/fe_lagrange_gll_shape_1D.C +31 +19 +38.71%
src/fe/fe_lagrange_gll_shape_2D.C +31 +19 +38.71%
src/fe/fe_lagrange_gll_shape_3D.C +31 +19 +38.71%
src/fe/fe_type.C +3 0 +0.99%
TOTAL +336 +111 +0.00%
code
coverage unchanged
code
coverage increased
code
coverage decreased
+
line added or modified

include/fe/fe_lagrange_gll_shape.h

52  
53  
54  
55 +
56  
57 +
58  
59  
60  

 * \returns The number of distinct second derivatives in dimension \p Dim.
 */
template <unsigned int Dim>
constexpr unsigned int fe_lagrange_gll_n_second_derivs()
{
  return Dim * (Dim + 1) / 2;
}

/**

81  
82  
83  
84 +
85  
86  
87  
88  
89  
90 +
91  
92 +
93  
94 +
95  
96 +
97 +
98  
99  
100 +
101 +
102  
103  
104  

 * element argument is unused.
 */
template <>
inline void fe_lagrange_gll_tensor_index<L2_LAGRANGE_GLL>(const Elem *,
                                                          const Order order,
                                                          const unsigned int dim,
                                                          const unsigned int i,
                                                          unsigned int * indices)
{
  const unsigned int n = static_cast<unsigned int>(order) + 1;

  unsigned int stride = 1;

  for (const auto d : make_range(dim))
    {
      indices[d] = (i / stride) % n;
      stride *= n;
    }

  libmesh_assert_less (i, stride);
}

/**
 * \returns Shape function \p i of family \p T in dimension \p Dim, differentiated once

110  
111  
112  
113 +
114  
115  
116  

 * FE::shape_deriv() passes one, and FE::shape_second_deriv() passes two.
 */
template <FEFamily T, unsigned int Dim>
inline Real fe_lagrange_gll_tensor_product(const Elem * elem,
                                           const Order order,
                                           const unsigned int i,
                                           const Point & p,

119  
120  
121  
122 +
123  
124 +
125  
126 +
127  
128 +
129 +
130 +
131 +
132  
133 +
134  
135 +
136 +
137 +
138  
139 +
140 +
141 +
142  
143  
144 +
145 +
146 +
147  
148  
149 +
150 +
151  
152  
153  
154 +
155  
156  
157  

{
  unsigned int indices[3];

  fe_lagrange_gll_tensor_index<T>(elem, order, Dim, i, indices);

  Real value = 1.;

  for (const auto d : make_range(Dim))
    {
      unsigned int deriv_order = 0;
      for (const auto k : make_range(n_deriv_coords))
        if (deriv_coords[k] == d)
          ++deriv_order;

      switch (deriv_order)
        {
        case 0:
          value *= fe_lagrange_gll_1D_shape(order, indices[d], p(d));
          break;

        case 1:
          value *= fe_lagrange_gll_1D_shape_deriv(order, indices[d], 0, p(d));
          break;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
        case 2:
          value *= fe_lagrange_gll_1D_shape_second_deriv(order, indices[d], 0, p(d));
          break;
#endif

        default:
          libmesh_error_msg("Invalid derivative order " << deriv_order << ".");
        }
    }

  return value;
}

} // namespace libMesh

include/fe/fe_lagrange_gll_shape_1D.h

78  
79  
80  
81 +
82  
83 +
84  
85  
86  
87 +
88 +
89 +
90  
91 +
92 +
93 +
94  
95 +
96 +
97  
98 +
99 +
100 +
101 +
102 +
103  
104  
105  
106  
107  
108  
109 +
110 +
111  
112 +
113 +
114 +
115  
116 +
117 +
118  
119 +
120  
121  
122  
123  
124  
125  
126 +
127 +
128  
129 +
130 +
131 +
132  
133 +
134 +
135 +
136  
137 +
138  
139  
140  
141 +
142 +
143  
144  
145  

/**
 * \returns The one-dimensional nodal basis on \p n Gauss-Lobatto points.
 */
inline GLLBasis1D build_gll_basis_1D(const unsigned int n)
{
  GLLBasis1D b;

  // QGaussLobatto chooses its number of points from the order: a rule of n points is exact
  // through degree 2n-3
  QGaussLobatto rule(1, static_cast<Order>(2*n - 3));
  rule.init(EDGE2);
  libmesh_assert_equal_to(rule.n_points(), n);

  b.nodes.resize(n);
  for (const auto i : make_range(n))
    b.nodes[i] = rule.qp(i)(0);

  b.weights.resize(n);
  for (const auto i : make_range(n))
    {
      Real product = 1.;
      for (const auto j : make_range(n))
        if (j != i)
          product *= b.nodes[i] - b.nodes[j];
      b.weights[i] = 1. / product;
    }

  // Entry (k,i) is the derivative of shape function i at point k. Off the diagonal it is
  // (w_i/w_k)/(x_k - x_i). The shape functions sum to one, so each row sums to zero, and
  // the diagonal entry is minus the sum of the rest of its row. That also gives the
  // endpoint entries -p(p+1)/4 and p(p+1)/4, so they don't need a special case.
  b.deriv.assign(n, std::vector<Real>(n, 0.));
  for (const auto k : make_range(n))
    {
      Real row_sum = 0.;
      for (const auto i : make_range(n))
        if (i != k)
          {
            b.deriv[k][i] = (b.weights[i] / b.weights[k]) / (b.nodes[k] - b.nodes[i]);
            row_sum += b.deriv[k][i];
          }
      b.deriv[k][k] = -row_sum;
    }

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
  // Off the diagonal, entry (k,i) of the second differentiation matrix is
  // 2 D_ki (D_kk - 1/(x_k - x_i)), with D the first differentiation matrix. Each row sums to
  // zero for the same reason as in D, which gives the diagonal.
  b.second.assign(n, std::vector<Real>(n, 0.));
  for (const auto k : make_range(n))
    {
      Real row_sum = 0.;
      for (const auto i : make_range(n))
        if (i != k)
          {
            b.second[k][i] =
              2. * b.deriv[k][i] * (b.deriv[k][k] - 1. / (b.nodes[k] - b.nodes[i]));
            row_sum += b.second[k][i];
          }
      b.second[k][k] = -row_sum;
    }
#endif

  return b;
}

/**
 * \returns The one-dimensional basis of degree \p order.

149  
150  
151  
152 +
153  
154 +
155  
156 +
157 +
158 +
159 +
160 +
161  
162 +
163  
164 +
165  
166  
167  
168  
169 +
170  
171  
172  

 * of n points costs \f$O(n^2)\f$ operations, and building all of them takes well under a
 * millisecond.
 */
inline const GLLBasis1D & fe_lagrange_gll_1D_basis(const Order order)
{
  static const std::array<GLLBasis1D, fe_lagrange_gll_max_points_1D + 1> bases = []()
    {
      std::array<GLLBasis1D, fe_lagrange_gll_max_points_1D + 1> b;
      for (const auto n : make_range(2u, fe_lagrange_gll_max_points_1D + 1))
        b[n] = build_gll_basis_1D(n);
      return b;
    }();

  const unsigned int n = static_cast<unsigned int>(order) + 1;

  libmesh_error_msg_if(n < 2 || n > fe_lagrange_gll_max_points_1D,
                       "The Gauss-Lobatto nodal basis supports orders 1 through "
                       << fe_lagrange_gll_max_points_1D - 1 << ", but order " << order
                       << " was requested.");

  return bases[n];
}



178  
179  
180  
181 +
182  
183 +
184  
185 +
186 +
187  
188 +
189  
190  
191  

 * is exact. The points of a collocated quadrature rule still match, because the basis and
 * the rule read the same values from QGaussLobatto.
 */
inline unsigned int fe_lagrange_gll_1D_point(const GLLBasis1D & b, const Real xi)
{
  const auto it = std::lower_bound(b.nodes.begin(), b.nodes.end(), xi);

  if (it != b.nodes.end() && *it == xi)
    return cast_int<unsigned int>(std::distance(b.nodes.begin(), it));

  return invalid_uint;
}

/**

196  
197  
198  
199 +
200  
201  
202  
203 +
204  
205 +
206  
207 +
208 +
209 +
210  
211 +
212  
213 +
214  
215 +
216  
217 +
218  
219 +
220 +
221  
222  
223 +
224  
225  
226  

 * relies on. Elsewhere it is evaluated with the barycentric formula
 * \f$\ell_i(\xi) = \frac{w_i/(\xi - \xi_i)}{\sum_j w_j/(\xi - \xi_j)}\f$.
 */
inline Real fe_lagrange_gll_1D_shape(const Order order,
                                     const unsigned int i,
                                     const Real xi)
{
  const GLLBasis1D & b = fe_lagrange_gll_1D_basis(order);

  libmesh_assert_less (i, b.nodes.size());

  const unsigned int k = fe_lagrange_gll_1D_point(b, xi);
  if (k != invalid_uint)
    return (i == k) ? 1. : 0.;

  Real numerator = 0., denominator = 0.;

  for (const auto j : index_range(b.nodes))
    {
      const Real term = b.weights[j] / (xi - b.nodes[j]);

      denominator += term;

      if (j == i)
        numerator = term;
    }

  return numerator / denominator;
}

/**

231  
232  
233  
234 +
235  
236  
237  
238 +
239  
240 +
241  
242 +
243  
244 +
245 +
246  
247  
248 +
249  
250  
251  

 * accurate as \p xi approaches a point.
 */
template <typename Column>
inline Real fe_lagrange_gll_1D_interpolate(const GLLBasis1D & b,
                                           const Real xi,
                                           const Column & f)
{
  Real numerator = 0., denominator = 0.;

  for (const auto k : index_range(b.nodes))
    {
      const Real term = b.weights[k] / (xi - b.nodes[k]);

      numerator += term * f(k);
      denominator += term;
    }

  return numerator / denominator;
}

/**

257  
258  
259  
260 +
261  
262  
263  
264  
265  
266 +
267  
268 +
269  
270 +
271  
272 +
273 +
274 +
275  
276 +
277 +
278  
279  
280  

 * interpolant of column \p i of that matrix. Differentiating the barycentric formula
 * directly would instead cancel catastrophically as \p xi approaches a point.
 */
inline Real fe_lagrange_gll_1D_shape_deriv(const Order order,
                                           const unsigned int i,
                                           const unsigned int libmesh_dbg_var(j),
                                           const Real xi)
{
  // only d()/dxi in 1D!
  libmesh_assert_equal_to (j, 0);

  const GLLBasis1D & b = fe_lagrange_gll_1D_basis(order);

  libmesh_assert_less (i, b.nodes.size());

  const unsigned int k = fe_lagrange_gll_1D_point(b, xi);
  if (k != invalid_uint)
    return b.deriv[k][i];

  return fe_lagrange_gll_1D_interpolate(b, xi, [&b, i](const std::size_t m)
                                        { return b.deriv[m][i]; });
}

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

285  
286  
287  
288 +
289  
290  
291  
292  
293  
294 +
295  
296 +
297  
298 +
299  
300 +
301 +
302 +
303  
304 +
305 +
306  
307  
308  

 * At a basis point it is read from the second differentiation matrix. Elsewhere it is the
 * interpolant of column \p i of that matrix, for the same reason as the first derivative.
 */
inline Real fe_lagrange_gll_1D_shape_second_deriv(const Order order,
                                                  const unsigned int i,
                                                  const unsigned int libmesh_dbg_var(j),
                                                  const Real xi)
{
  // only d^2()/dxi^2 in 1D!
  libmesh_assert_equal_to (j, 0);

  const GLLBasis1D & b = fe_lagrange_gll_1D_basis(order);

  libmesh_assert_less (i, b.nodes.size());

  const unsigned int k = fe_lagrange_gll_1D_point(b, xi);
  if (k != invalid_uint)
    return b.second[k][i];

  return fe_lagrange_gll_1D_interpolate(b, xi, [&b, i](const std::size_t m)
                                        { return b.second[m][i]; });
}
#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

include/fe/fe_lagrange_shape_1D.h

250  
251  
252  
253  
254  
255  
256  
257  
258  
259  
260  
261  
262  
263  
264  
265  
266  
267  
268  
269  
270  
271  
272  
273  
274  
275  
276  
277  
278  



inline
Real fe_lagrange_1D_cubic_shape_second_deriv(const unsigned int i,
                                             const unsigned int libmesh_dbg_var(j),
                                             const Real xi)
{
  // Don't need to switch on j.  1D shape functions
  // depend on xi only!
  libmesh_assert_equal_to (j, 0);
  libmesh_assert_less(i, 4);

  switch (i)
    {
    case 0:
      return -9./16.*(6.*xi-2);

    case 1:
      return -9./16.*(-6*xi-2.);

    case 2:
      return 27./16.*(6*xi-2./3.);

    // case 2
    default:
      return 27./16.*(-6*xi-2./3.);
    }
}


296  
297  
298  
299  
300  
301  
302  
303  

      return fe_lagrange_1D_quadratic_shape_second_deriv(i, j, xi);

    // case THIRD
    default:
      return fe_lagrange_1D_cubic_shape_second_deriv(i, j, xi);
    } // end switch (order)
}

src/fe/fe_abstract.C

100  
101  
102  
103 +
104 +
105  
106  
107  

          case L2_LAGRANGE:
            return std::make_unique<FE<0,L2_LAGRANGE>>(fet);

          case L2_LAGRANGE_GLL:
            return std::make_unique<FE<0,L2_LAGRANGE_GLL>>(fet);

          case L2_LAGRANGE_VEC:
            return std::make_unique<FE<0,L2_LAGRANGE_VEC>>(fet);

180  
181  
182  
183 +
184 +
185  
186  
187  

          case L2_LAGRANGE:
            return std::make_unique<FE<1,L2_LAGRANGE>>(fet);

          case L2_LAGRANGE_GLL:
            return std::make_unique<FE<1,L2_LAGRANGE_GLL>>(fet);

          case L2_LAGRANGE_VEC:
            return std::make_unique<FE<1,L2_LAGRANGE_VEC>>(fet);

262  
263  
264  
265 +
266 +
267  
268  
269  

          case L2_LAGRANGE:
            return std::make_unique<FE<2,L2_LAGRANGE>>(fet);

          case L2_LAGRANGE_GLL:
            return std::make_unique<FE<2,L2_LAGRANGE_GLL>>(fet);

          case L2_LAGRANGE_VEC:
            return std::make_unique<FE<2,L2_LAGRANGE_VEC>>(fet);

344  
345  
346  
347 +
348 +
349  
350  
351  

          case L2_LAGRANGE:
            return std::make_unique<FE<3,L2_LAGRANGE>>(fet);

          case L2_LAGRANGE_GLL:
            return std::make_unique<FE<3,L2_LAGRANGE_GLL>>(fet);

          case L2_LAGRANGE_VEC:
            return std::make_unique<FE<3,L2_LAGRANGE_VEC>>(fet);

435  
436  
437  
438  
439  
440  
441  
442  
443 +
444  
445  
446  
447  

        nodes[1] = Point (1.,0.,0.);
        return;
      }
    case EDGE4: // not nested with EDGE3
      {
        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.);
        return;
      }
    case TRI7:
      {

src/fe/fe_base.C

210  
211  
212  
213 +
214 +
215  
216  
217  

          case L2_LAGRANGE:
            return std::make_unique<FE<0,L2_LAGRANGE>>(fet);

          case L2_LAGRANGE_GLL:
            return std::make_unique<FE<0,L2_LAGRANGE_GLL>>(fet);

          case HIERARCHIC:
            return std::make_unique<FE<0,HIERARCHIC>>(fet);

263  
264  
265  
266 +
267 +
268  
269  
270  

          case L2_LAGRANGE:
            return std::make_unique<FE<1,L2_LAGRANGE>>(fet);

          case L2_LAGRANGE_GLL:
            return std::make_unique<FE<1,L2_LAGRANGE_GLL>>(fet);

          case HIERARCHIC:
            return std::make_unique<FE<1,HIERARCHIC>>(fet);

318  
319  
320  
321 +
322 +
323  
324  
325  

          case L2_LAGRANGE:
            return std::make_unique<FE<2,L2_LAGRANGE>>(fet);

          case L2_LAGRANGE_GLL:
            return std::make_unique<FE<2,L2_LAGRANGE_GLL>>(fet);

          case HIERARCHIC:
            return std::make_unique<FE<2,HIERARCHIC>>(fet);

376  
377  
378  
379 +
380 +
381  
382  
383  

          case L2_LAGRANGE:
            return std::make_unique<FE<3,L2_LAGRANGE>>(fet);

          case L2_LAGRANGE_GLL:
            return std::make_unique<FE<3,L2_LAGRANGE_GLL>>(fet);

          case HIERARCHIC:
            return std::make_unique<FE<3,HIERARCHIC>>(fet);

src/fe/fe_boundary.C

79  
80  
81  
82 +
83  
84  
85  

LIBMESH_ERRORS_IN_LOW_D(SIDE_HIERARCHIC)
LIBMESH_ERRORS_IN_LOW_D(LAGRANGE)
LIBMESH_ERRORS_IN_LOW_D(L2_LAGRANGE)
LIBMESH_ERRORS_IN_LOW_D(L2_LAGRANGE_GLL)
LIBMESH_ERRORS_IN_LOW_D(LAGRANGE_VEC)
LIBMESH_ERRORS_IN_LOW_D(L2_LAGRANGE_VEC)
LIBMESH_ERRORS_IN_LOW_D(MONOMIAL)

src/fe/fe_interface.C

2203  
2204  
2205  
2206 +
2207 +
2208  
2209 +
2210  
2211  
2212  

          return unknown;
        }
      break;
    case L2_LAGRANGE_GLL:
      switch (el_t)
        {
        case EDGE2:
        case EDGE3:
        case EDGE4:
        case QUAD4:

2220  
2221  
2222  
2223 +
2224 +
2225 +
2226  
2227  
2228  

        case HEX27:
          // A basis of degree p needs p+1 Gauss-Lobatto points, and QGaussLobatto tabulates
          // at most fe_lagrange_gll_max_points_1D
          return fe_lagrange_gll_max_points_1D - 1;
        default:
          return unknown;
        }
      break;
    case MONOMIAL:

src/fe/fe_l2_lagrange_gll.C

38  
39  
40  
41 +
42  
43 +
44  
45 +
46  
47 +
48  
49  
50  
51  
52 +
53  
54 +
55 +
56  
57 +
58  
59  
60 +
61  
62 +
63  
64  
65  
66  
67  
68 +
69  
70 +
71  
72  
73 +
74  
75 +
76 +
77  
78 +
79 +
80  
81  
82  
83  
84  
85 +
86  
87 +
88 +
89  
90  
91  

 * to the element, so the count doesn't depend on the number of nodes, and every supported
 * order works on a first order element.
 */
unsigned int l2_lagrange_gll_n_dofs(const ElemType t, const Order o)
{
  libmesh_assert_greater (o, 0);

  const unsigned int n = static_cast<unsigned int>(o) + 1;

  libmesh_error_msg_if(n > fe_lagrange_gll_max_points_1D,
                       "The L2_LAGRANGE_GLL basis of order " << o << " needs " << n
                       << " Gauss-Lobatto points, but QGaussLobatto tabulates at most "
                       << fe_lagrange_gll_max_points_1D << ".");

  switch (t)
    {
    case NODEELEM:
      return 1;

    case EDGE2:
    case EDGE3:
    case EDGE4:
      return n;

    case QUAD4:
    case QUADSHELL4:
    case QUAD8:
    case QUADSHELL8:
    case QUAD9:
    case QUADSHELL9:
      return n*n;

    case HEX8:
    case HEX20:
    case HEX27:
      return n*n*n;

    case INVALID_ELEM:
      return 0;

    default:
      libmesh_error_msg("The L2_LAGRANGE_GLL basis is a tensor product and is not defined on "
                        << Utility::enum_to_string(t) << ".");
    }
}


unsigned int l2_lagrange_gll_n_dofs(const Elem * e, const Order o)
{
  libmesh_assert(e);
  return l2_lagrange_gll_n_dofs(e->type(), o);
}

/**

95  
96  
97  
98 +
99  
100  
101  
102  
103  
104 +
105  
106 +
107  
108 +
109  
110  
111 +
112  
113 +
114 +
115 +
116 +
117  
118 +
119  
120 +
121 +
122 +
123 +
124 +
125  
126 +
127 +
128  
129  
130  
131 +
132 +
133 +
134 +
135  
136 +
137 +
138 +
139 +
140  
141 +
142 +
143 +
144 +
145  
146 +
147 +
148 +
149 +
150  
151 +
152 +
153 +
154 +
155  
156 +
157 +
158 +
159 +
160  
161  
162 +
163 +
164 +
165 +
166  
167  
168 +
169 +
170 +
171 +
172  
173  
174 +
175 +
176 +
177 +
178  
179  
180  
181  
182 +
183  
184  
185  
186 +
187  
188  
189 +
190  
191  
192  
193 +
194  
195  
196  

 * The interpolation points coincide with the element nodes only at orders one and two, so
 * the solution is evaluated through the shape functions at every order.
 */
void l2_lagrange_gll_nodal_soln(const Elem * elem,
                                const Order order,
                                const std::vector<Number> & elem_soln,
                                std::vector<Number> & nodal_soln,
                                const bool add_p_level)
{
  const unsigned int n_nodes = elem->n_nodes();

  nodal_soln.resize(n_nodes);

  const FEType fe_type(order, L2_LAGRANGE_GLL);

  const unsigned int n_sf =
    FEInterface::n_shape_functions(fe_type, elem, add_p_level);

  std::vector<Point> refspace_nodes;
  FEBase::get_refspace_nodes(elem->type(), refspace_nodes);
  libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
  libmesh_assert_equal_to (elem_soln.size(), n_sf);

  std::fill(nodal_soln.begin(), nodal_soln.end(), 0);

  for (const auto n : make_range(n_nodes))
    for (const auto i : make_range(n_sf))
      nodal_soln[n] += elem_soln[i] *
        FEInterface::shape(fe_type, elem, i, refspace_nodes[n], add_p_level);
}

LIBMESH_FE_NODAL_SOLN(L2_LAGRANGE_GLL, l2_lagrange_gll_nodal_soln)
LIBMESH_FE_SIDE_NODAL_SOLN(L2_LAGRANGE_GLL)


// Every degree of freedom of this basis belongs to the element
template <> unsigned int FE<0,L2_LAGRANGE_GLL>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_gll_n_dofs(t, o); }
template <> unsigned int FE<1,L2_LAGRANGE_GLL>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_gll_n_dofs(t, o); }
template <> unsigned int FE<2,L2_LAGRANGE_GLL>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_gll_n_dofs(t, o); }
template <> unsigned int FE<3,L2_LAGRANGE_GLL>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_gll_n_dofs(t, o); }

template <> unsigned int FE<0,L2_LAGRANGE_GLL>::n_dofs(const Elem * e, const Order o) { return l2_lagrange_gll_n_dofs(e, o); }
template <> unsigned int FE<1,L2_LAGRANGE_GLL>::n_dofs(const Elem * e, const Order o) { return l2_lagrange_gll_n_dofs(e, o); }
template <> unsigned int FE<2,L2_LAGRANGE_GLL>::n_dofs(const Elem * e, const Order o) { return l2_lagrange_gll_n_dofs(e, o); }
template <> unsigned int FE<3,L2_LAGRANGE_GLL>::n_dofs(const Elem * e, const Order o) { return l2_lagrange_gll_n_dofs(e, o); }

template <> unsigned int FE<0,L2_LAGRANGE_GLL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
template <> unsigned int FE<1,L2_LAGRANGE_GLL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
template <> unsigned int FE<2,L2_LAGRANGE_GLL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
template <> unsigned int FE<3,L2_LAGRANGE_GLL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }

template <> unsigned int FE<0,L2_LAGRANGE_GLL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
template <> unsigned int FE<1,L2_LAGRANGE_GLL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
template <> unsigned int FE<2,L2_LAGRANGE_GLL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
template <> unsigned int FE<3,L2_LAGRANGE_GLL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }

template <> unsigned int FE<0,L2_LAGRANGE_GLL>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_gll_n_dofs(t, o); }
template <> unsigned int FE<1,L2_LAGRANGE_GLL>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_gll_n_dofs(t, o); }
template <> unsigned int FE<2,L2_LAGRANGE_GLL>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_gll_n_dofs(t, o); }
template <> unsigned int FE<3,L2_LAGRANGE_GLL>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_gll_n_dofs(t, o); }

template <> unsigned int FE<0,L2_LAGRANGE_GLL>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_lagrange_gll_n_dofs(&e, o); }
template <> unsigned int FE<1,L2_LAGRANGE_GLL>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_lagrange_gll_n_dofs(&e, o); }
template <> unsigned int FE<2,L2_LAGRANGE_GLL>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_lagrange_gll_n_dofs(&e, o); }
template <> unsigned int FE<3,L2_LAGRANGE_GLL>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_lagrange_gll_n_dofs(&e, o); }


template <> FEContinuity FE<0,L2_LAGRANGE_GLL>::get_continuity() const { return DISCONTINUOUS; }
template <> FEContinuity FE<1,L2_LAGRANGE_GLL>::get_continuity() const { return DISCONTINUOUS; }
template <> FEContinuity FE<2,L2_LAGRANGE_GLL>::get_continuity() const { return DISCONTINUOUS; }
template <> FEContinuity FE<3,L2_LAGRANGE_GLL>::get_continuity() const { return DISCONTINUOUS; }

// A nodal basis is not hierarchic
template <> bool FE<0,L2_LAGRANGE_GLL>::is_hierarchic() const { return false; }
template <> bool FE<1,L2_LAGRANGE_GLL>::is_hierarchic() const { return false; }
template <> bool FE<2,L2_LAGRANGE_GLL>::is_hierarchic() const { return false; }
template <> bool FE<3,L2_LAGRANGE_GLL>::is_hierarchic() const { return false; }

// The shape functions are the same on every element of a type
template <> bool FE<0,L2_LAGRANGE_GLL>::shapes_need_reinit() const { return false; }
template <> bool FE<1,L2_LAGRANGE_GLL>::shapes_need_reinit() const { return false; }
template <> bool FE<2,L2_LAGRANGE_GLL>::shapes_need_reinit() const { return false; }
template <> bool FE<3,L2_LAGRANGE_GLL>::shapes_need_reinit() const { return false; }

// This DISCONTINUOUS basis doesn't need hanging node constraints, so these are empty
#ifdef LIBMESH_ENABLE_AMR
template <>
void FE<2,L2_LAGRANGE_GLL>::compute_constraints (DofConstraints &,
                                                 DofMap &,
                                                 const unsigned int,
                                                 const Elem *)
{ }

template <>
void FE<3,L2_LAGRANGE_GLL>::compute_constraints (DofConstraints &,
                                                 DofMap &,
                                                 const unsigned int,
                                                 const Elem *)
{ }
#endif // LIBMESH_ENABLE_AMR

src/fe/fe_lagrange_gll_shape_0D.C

26  
27  
28  
29 +
30  
31  
32  
33 +
34  
35  
36  
37  
38 +
39 +
40  
41  
42  
43  
44 +
45  
46  
47  
48  
49  
50 +
51 +
52  
53  
54  
55  
56 +
57  
58  
59  
60  
61  
62 +
63 +
64  
65  
66  
67  
68 +
69  
70  
71  
72  
73  
74 +
75  
76  
77  
78  
79  
80 +
81  
82  
83  
84  
85  
86  
87 +
88  
89  
90  
91  
92  
93 +
94  
95  
96  
97  
98  
99  
100 +
101  
102  
103  

{


LIBMESH_DEFAULT_VECTORIZED_FE(0,L2_LAGRANGE_GLL)


template <>
Real FE<0,L2_LAGRANGE_GLL>::shape(const ElemType,
                                  const Order,
                                  const unsigned int libmesh_dbg_var(i),
                                  const Point &)
{
  libmesh_assert_less (i, 1);
  return 1.;
}


template <>
Real FE<0,L2_LAGRANGE_GLL>::shape(const Elem *,
                                  const Order,
                                  const unsigned int libmesh_dbg_var(i),
                                  const Point &,
                                  const bool)
{
  libmesh_assert_less (i, 1);
  return 1.;
}


template <>
Real FE<0,L2_LAGRANGE_GLL>::shape(const FEType,
                                  const Elem *,
                                  const unsigned int libmesh_dbg_var(i),
                                  const Point &,
                                  const bool)
{
  libmesh_assert_less (i, 1);
  return 1.;
}


template <>
Real FE<0,L2_LAGRANGE_GLL>::shape_deriv(const ElemType,
                                        const Order,
                                        const unsigned int,
                                        const unsigned int,
                                        const Point &)
{
  libmesh_error_msg("No spatial derivatives in 0D!");
  return 0.;
}


template <>
Real FE<0,L2_LAGRANGE_GLL>::shape_deriv(const Elem *,
                                        const Order,
                                        const unsigned int,
                                        const unsigned int,
                                        const Point &,
                                        const bool)
{
  libmesh_error_msg("No spatial derivatives in 0D!");
  return 0.;
}


template <>
Real FE<0,L2_LAGRANGE_GLL>::shape_deriv(const FEType,
                                        const Elem *,
                                        const unsigned int,
                                        const unsigned int,
                                        const Point &,
                                        const bool)
{
  libmesh_error_msg("No spatial derivatives in 0D!");
  return 0.;
}


105  
106  
107  
108 +
109  
110  
111  
112  
113  
114 +
115  
116  
117  
118  
119  
120 +
121  
122  
123  
124  
125  
126  
127 +
128  
129  
130  
131  
132  
133 +
134  
135  
136  
137  
138  
139  
140 +
141  
142  
143  

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

template <>
Real FE<0,L2_LAGRANGE_GLL>::shape_second_deriv(const ElemType,
                                               const Order,
                                               const unsigned int,
                                               const unsigned int,
                                               const Point &)
{
  libmesh_error_msg("No spatial derivatives in 0D!");
  return 0.;
}


template <>
Real FE<0,L2_LAGRANGE_GLL>::shape_second_deriv(const Elem *,
                                               const Order,
                                               const unsigned int,
                                               const unsigned int,
                                               const Point &,
                                               const bool)
{
  libmesh_error_msg("No spatial derivatives in 0D!");
  return 0.;
}


template <>
Real FE<0,L2_LAGRANGE_GLL>::shape_second_deriv(const FEType,
                                               const Elem *,
                                               const unsigned int,
                                               const unsigned int,
                                               const Point &,
                                               const bool)
{
  libmesh_error_msg("No spatial derivatives in 0D!");
  return 0.;
}

src/fe/fe_lagrange_gll_shape_1D.C

27  
28  
29  
30 +
31  
32  
33  
34 +
35  
36  
37  
38  
39  
40 +
41  
42  
43  
44  
45 +
46  
47  
48  
49  
50  
51 +
52  
53  
54 +
55  
56  
57  
58  
59 +
60  
61  
62  
63  
64  
65 +
66  
67  
68 +
69  
70  
71  
72  
73 +
74  
75  
76  
77  
78  
79 +
80  
81  
82 +
83  
84  
85  
86  
87 +
88  
89  
90  
91  
92  
93  
94 +
95 +
96  
97  
98 +
99  
100  
101  
102  
103 +
104  
105  
106  
107  
108  
109  
110 +
111 +
112  
113  
114 +
115  
116  
117  
118  
119  
120  
121 +
122  
123  
124  
125  
126  
127 +
128  
129  
130 +
131  
132  
133  
134  
135 +
136  
137  
138  
139  
140  
141  
142 +
143 +
144  
145  
146 +
147  
148  
149  
150  
151 +
152  
153  
154  
155  
156  
157  
158 +
159 +
160  
161  
162 +
163  
164  
165  

{


LIBMESH_DEFAULT_VECTORIZED_FE(1,L2_LAGRANGE_GLL)


template <>
Real FE<1,L2_LAGRANGE_GLL>::shape(const ElemType,
                                  const Order order,
                                  const unsigned int i,
                                  const Point & p)
{
  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1>
    (nullptr, order, i, p, nullptr, 0);
}


template <>
Real FE<1,L2_LAGRANGE_GLL>::shape(const Elem * elem,
                                  const Order order,
                                  const unsigned int i,
                                  const Point & p,
                                  const bool add_p_level)
{
  libmesh_assert(elem);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1>
    (elem, order + add_p_level*elem->p_level(), i, p, nullptr, 0);
}


template <>
Real FE<1,L2_LAGRANGE_GLL>::shape(const FEType fet,
                                  const Elem * elem,
                                  const unsigned int i,
                                  const Point & p,
                                  const bool add_p_level)
{
  libmesh_assert(elem);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1>
    (elem, fet.order + add_p_level*elem->p_level(), i, p, nullptr, 0);
}


template <>
Real FE<1,L2_LAGRANGE_GLL>::shape_deriv(const ElemType,
                                        const Order order,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p)
{
  libmesh_assert_less (j, 1);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1>
    (nullptr, order, i, p, &j, 1);
}


template <>
Real FE<1,L2_LAGRANGE_GLL>::shape_deriv(const Elem * elem,
                                        const Order order,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p,
                                        const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, 1);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1>
    (elem, order + add_p_level*elem->p_level(), i, p, &j, 1);
}


template <>
Real FE<1,L2_LAGRANGE_GLL>::shape_deriv(const FEType fet,
                                        const Elem * elem,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p,
                                        const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, 1);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1>
    (elem, fet.order + add_p_level*elem->p_level(), i, p, &j, 1);
}


#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

template <>
Real FE<1,L2_LAGRANGE_GLL>::shape_second_deriv(const ElemType,
                                               const Order order,
                                               const unsigned int i,
                                               const unsigned int j,
                                               const Point & p)
{
  libmesh_assert_less (j, fe_lagrange_gll_n_second_derivs<1>());

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1>
    (nullptr, order, i, p, fe_lagrange_gll_second_deriv_pairs[j], 2);
}


template <>
Real FE<1,L2_LAGRANGE_GLL>::shape_second_deriv(const Elem * elem,
                                               const Order order,
                                               const unsigned int i,
                                               const unsigned int j,
                                               const Point & p,
                                               const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, fe_lagrange_gll_n_second_derivs<1>());

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1>
    (elem, order + add_p_level*elem->p_level(), i, p, fe_lagrange_gll_second_deriv_pairs[j], 2);
}


template <>
Real FE<1,L2_LAGRANGE_GLL>::shape_second_deriv(const FEType fet,
                                               const Elem * elem,
                                               const unsigned int i,
                                               const unsigned int j,
                                               const Point & p,
                                               const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, fe_lagrange_gll_n_second_derivs<1>());

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 1>
    (elem, fet.order + add_p_level*elem->p_level(), i, p, fe_lagrange_gll_second_deriv_pairs[j], 2);
}

src/fe/fe_lagrange_gll_shape_2D.C

27  
28  
29  
30 +
31  
32  
33  
34 +
35  
36  
37  
38  
39  
40 +
41  
42  
43  
44  
45 +
46  
47  
48  
49  
50  
51 +
52  
53  
54 +
55  
56  
57  
58  
59 +
60  
61  
62  
63  
64  
65 +
66  
67  
68 +
69  
70  
71  
72  
73 +
74  
75  
76  
77  
78  
79 +
80  
81  
82 +
83  
84  
85  
86  
87 +
88  
89  
90  
91  
92  
93  
94 +
95 +
96  
97  
98 +
99  
100  
101  
102  
103 +
104  
105  
106  
107  
108  
109  
110 +
111 +
112  
113  
114 +
115  
116  
117  
118  
119  
120  
121 +
122  
123  
124  
125  
126  
127 +
128  
129  
130 +
131  
132  
133  
134  
135 +
136  
137  
138  
139  
140  
141  
142 +
143 +
144  
145  
146 +
147  
148  
149  
150  
151 +
152  
153  
154  
155  
156  
157  
158 +
159 +
160  
161  
162 +
163  
164  
165  

{


LIBMESH_DEFAULT_VECTORIZED_FE(2,L2_LAGRANGE_GLL)


template <>
Real FE<2,L2_LAGRANGE_GLL>::shape(const ElemType,
                                  const Order order,
                                  const unsigned int i,
                                  const Point & p)
{
  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 2>
    (nullptr, order, i, p, nullptr, 0);
}


template <>
Real FE<2,L2_LAGRANGE_GLL>::shape(const Elem * elem,
                                  const Order order,
                                  const unsigned int i,
                                  const Point & p,
                                  const bool add_p_level)
{
  libmesh_assert(elem);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 2>
    (elem, order + add_p_level*elem->p_level(), i, p, nullptr, 0);
}


template <>
Real FE<2,L2_LAGRANGE_GLL>::shape(const FEType fet,
                                  const Elem * elem,
                                  const unsigned int i,
                                  const Point & p,
                                  const bool add_p_level)
{
  libmesh_assert(elem);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 2>
    (elem, fet.order + add_p_level*elem->p_level(), i, p, nullptr, 0);
}


template <>
Real FE<2,L2_LAGRANGE_GLL>::shape_deriv(const ElemType,
                                        const Order order,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p)
{
  libmesh_assert_less (j, 2);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 2>
    (nullptr, order, i, p, &j, 1);
}


template <>
Real FE<2,L2_LAGRANGE_GLL>::shape_deriv(const Elem * elem,
                                        const Order order,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p,
                                        const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, 2);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 2>
    (elem, order + add_p_level*elem->p_level(), i, p, &j, 1);
}


template <>
Real FE<2,L2_LAGRANGE_GLL>::shape_deriv(const FEType fet,
                                        const Elem * elem,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p,
                                        const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, 2);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 2>
    (elem, fet.order + add_p_level*elem->p_level(), i, p, &j, 1);
}


#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

template <>
Real FE<2,L2_LAGRANGE_GLL>::shape_second_deriv(const ElemType,
                                               const Order order,
                                               const unsigned int i,
                                               const unsigned int j,
                                               const Point & p)
{
  libmesh_assert_less (j, fe_lagrange_gll_n_second_derivs<2>());

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 2>
    (nullptr, order, i, p, fe_lagrange_gll_second_deriv_pairs[j], 2);
}


template <>
Real FE<2,L2_LAGRANGE_GLL>::shape_second_deriv(const Elem * elem,
                                               const Order order,
                                               const unsigned int i,
                                               const unsigned int j,
                                               const Point & p,
                                               const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, fe_lagrange_gll_n_second_derivs<2>());

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 2>
    (elem, order + add_p_level*elem->p_level(), i, p, fe_lagrange_gll_second_deriv_pairs[j], 2);
}


template <>
Real FE<2,L2_LAGRANGE_GLL>::shape_second_deriv(const FEType fet,
                                               const Elem * elem,
                                               const unsigned int i,
                                               const unsigned int j,
                                               const Point & p,
                                               const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, fe_lagrange_gll_n_second_derivs<2>());

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 2>
    (elem, fet.order + add_p_level*elem->p_level(), i, p, fe_lagrange_gll_second_deriv_pairs[j], 2);
}

src/fe/fe_lagrange_gll_shape_3D.C

27  
28  
29  
30 +
31  
32  
33  
34 +
35  
36  
37  
38  
39  
40 +
41  
42  
43  
44  
45 +
46  
47  
48  
49  
50  
51 +
52  
53  
54 +
55  
56  
57  
58  
59 +
60  
61  
62  
63  
64  
65 +
66  
67  
68 +
69  
70  
71  
72  
73 +
74  
75  
76  
77  
78  
79 +
80  
81  
82 +
83  
84  
85  
86  
87 +
88  
89  
90  
91  
92  
93  
94 +
95 +
96  
97  
98 +
99  
100  
101  
102  
103 +
104  
105  
106  
107  
108  
109  
110 +
111 +
112  
113  
114 +
115  
116  
117  
118  
119  
120  
121 +
122  
123  
124  
125  
126  
127 +
128  
129  
130 +
131  
132  
133  
134  
135 +
136  
137  
138  
139  
140  
141  
142 +
143 +
144  
145  
146 +
147  
148  
149  
150  
151 +
152  
153  
154  
155  
156  
157  
158 +
159 +
160  
161  
162 +
163  
164  
165  

{


LIBMESH_DEFAULT_VECTORIZED_FE(3,L2_LAGRANGE_GLL)


template <>
Real FE<3,L2_LAGRANGE_GLL>::shape(const ElemType,
                                  const Order order,
                                  const unsigned int i,
                                  const Point & p)
{
  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 3>
    (nullptr, order, i, p, nullptr, 0);
}


template <>
Real FE<3,L2_LAGRANGE_GLL>::shape(const Elem * elem,
                                  const Order order,
                                  const unsigned int i,
                                  const Point & p,
                                  const bool add_p_level)
{
  libmesh_assert(elem);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 3>
    (elem, order + add_p_level*elem->p_level(), i, p, nullptr, 0);
}


template <>
Real FE<3,L2_LAGRANGE_GLL>::shape(const FEType fet,
                                  const Elem * elem,
                                  const unsigned int i,
                                  const Point & p,
                                  const bool add_p_level)
{
  libmesh_assert(elem);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 3>
    (elem, fet.order + add_p_level*elem->p_level(), i, p, nullptr, 0);
}


template <>
Real FE<3,L2_LAGRANGE_GLL>::shape_deriv(const ElemType,
                                        const Order order,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p)
{
  libmesh_assert_less (j, 3);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 3>
    (nullptr, order, i, p, &j, 1);
}


template <>
Real FE<3,L2_LAGRANGE_GLL>::shape_deriv(const Elem * elem,
                                        const Order order,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p,
                                        const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, 3);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 3>
    (elem, order + add_p_level*elem->p_level(), i, p, &j, 1);
}


template <>
Real FE<3,L2_LAGRANGE_GLL>::shape_deriv(const FEType fet,
                                        const Elem * elem,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p,
                                        const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, 3);

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 3>
    (elem, fet.order + add_p_level*elem->p_level(), i, p, &j, 1);
}


#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

template <>
Real FE<3,L2_LAGRANGE_GLL>::shape_second_deriv(const ElemType,
                                               const Order order,
                                               const unsigned int i,
                                               const unsigned int j,
                                               const Point & p)
{
  libmesh_assert_less (j, fe_lagrange_gll_n_second_derivs<3>());

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 3>
    (nullptr, order, i, p, fe_lagrange_gll_second_deriv_pairs[j], 2);
}


template <>
Real FE<3,L2_LAGRANGE_GLL>::shape_second_deriv(const Elem * elem,
                                               const Order order,
                                               const unsigned int i,
                                               const unsigned int j,
                                               const Point & p,
                                               const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, fe_lagrange_gll_n_second_derivs<3>());

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 3>
    (elem, order + add_p_level*elem->p_level(), i, p, fe_lagrange_gll_second_deriv_pairs[j], 2);
}


template <>
Real FE<3,L2_LAGRANGE_GLL>::shape_second_deriv(const FEType fet,
                                               const Elem * elem,
                                               const unsigned int i,
                                               const unsigned int j,
                                               const Point & p,
                                               const bool add_p_level)
{
  libmesh_assert(elem);
  libmesh_assert_less (j, fe_lagrange_gll_n_second_derivs<3>());

  return fe_lagrange_gll_tensor_product<L2_LAGRANGE_GLL, 3>
    (elem, fet.order + add_p_level*elem->p_level(), i, p, fe_lagrange_gll_second_deriv_pairs[j], 2);
}

src/fe/fe_type.C

56  
57  
58  
59 +
60  
61 +
62 +
63  
64  
65  

  // p refinement preserves the pairing. QBase raises the order of a rule by twice the
  // element's p_level, and 2p-1 also rises by two per degree, so on an element of p level l
  // the rule has the p+l+1 points of the degree p+l basis.
  if (family == L2_LAGRANGE_GLL)
    {
      const int p = static_cast<int>(order.get_order());
      return std::make_unique<QGaussLobatto>(dim, static_cast<Order>(2*p - 1 + extraorder));
    }

  return std::make_unique<QGauss>(dim, static_cast<Order>(this->default_quadrature_order() + extraorder));