fe_lagrange_gll_shape.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
#ifndef LIBMESH_FE_LAGRANGE_GLL_SHAPE_H
#define LIBMESH_FE_LAGRANGE_GLL_SHAPE_H
 
// Local includes
#include "libmesh/enum_fe_family.h"
#include "libmesh/enum_order.h"
#include "libmesh/fe_lagrange_gll_shape_1D.h"
#include "libmesh/int_range.h"
 
// Tensor product evaluation shared by the Gauss-Lobatto nodal families. The families differ
// only in how a shape function index maps to a tensor index, so each family specializes
// fe_lagrange_gll_tensor_index and shares the rest.
 
namespace libMesh
{
 
constexpr unsigned int fe_lagrange_gll_second_deriv_pairs[6][2] =
  {{0, 0},   // d^2 phi / dxi^2
   {0, 1},   // d^2 phi / dxi deta
   {1, 1},   // d^2 phi / deta^2
   {0, 2},   // d^2 phi / dxi dzeta
   {1, 2},   // d^2 phi / deta dzeta
   {2, 2}};  // d^2 phi / dzeta^2
 
 
template <unsigned int Dim>
constexpr unsigned int fe_lagrange_gll_n_second_derivs()
{
  return Dim * (Dim + 1) / 2;
}
 
template <FEFamily T>
void fe_lagrange_gll_tensor_index(const Elem * elem,
                                  const Order order,
                                  const unsigned int dim,
                                  const unsigned int i,
                                  unsigned int * indices);
 
 
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);
}
 
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,
                                           const unsigned int * const deriv_coords,
                                           const unsigned int n_deriv_coords)
{
  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
 
#endif // LIBMESH_FE_LAGRANGE_GLL_SHAPE_H