fe_values_base.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_fe_values_base_h
#define dealii_fe_values_base_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/std_cxx20/iota_view.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/symmetric_tensor.h>
 
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/fe/fe_values_views.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/mapping_related_data.h>
 
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
 
#include <deal.II/hp/q_collection.h>
 
#include <deal.II/lac/read_vector.h>
 
#include <algorithm>
#include <memory>
#include <optional>
#include <type_traits>
#include <variant>
 
DEAL_II_NAMESPACE_OPEN
 
template <int dim, int spacedim>
class FEValuesBase : public Subscriptor
{
public:
  static constexpr unsigned int dimension = dim;
 
  static constexpr unsigned int space_dimension = spacedim;
 
  const unsigned int n_quadrature_points;
 
  const unsigned int max_n_quadrature_points;
 
  const unsigned int dofs_per_cell;
 
 
  FEValuesBase(const unsigned int                  n_q_points,
               const unsigned int                  dofs_per_cell,
               const UpdateFlags                   update_flags,
               const Mapping<dim, spacedim>       &mapping,
               const FiniteElement<dim, spacedim> &fe);
 
  FEValuesBase &
  operator=(const FEValuesBase &) = delete;
 
  FEValuesBase(const FEValuesBase &) = delete;
 
  virtual ~FEValuesBase() override;
 
  void
  always_allow_check_for_cell_similarity(const bool allow);
 
 
  const double &
  shape_value(const unsigned int i, const unsigned int q_point) const;
 
  double
  shape_value_component(const unsigned int i,
                        const unsigned int q_point,
                        const unsigned int component) const;
 
  const Tensor<1, spacedim> &
  shape_grad(const unsigned int i, const unsigned int q_point) const;
 
  Tensor<1, spacedim>
  shape_grad_component(const unsigned int i,
                       const unsigned int q_point,
                       const unsigned int component) const;
 
  const Tensor<2, spacedim> &
  shape_hessian(const unsigned int i, const unsigned int q_point) const;
 
  Tensor<2, spacedim>
  shape_hessian_component(const unsigned int i,
                          const unsigned int q_point,
                          const unsigned int component) const;
 
  const Tensor<3, spacedim> &
  shape_3rd_derivative(const unsigned int i, const unsigned int q_point) const;
 
  Tensor<3, spacedim>
  shape_3rd_derivative_component(const unsigned int i,
                                 const unsigned int q_point,
                                 const unsigned int component) const;
 
 
  template <typename Number>
  void
  get_function_values(const ReadVector<Number> &fe_function,
                      std::vector<Number>      &values) const;
 
  template <typename Number>
  void
  get_function_values(const ReadVector<Number>    &fe_function,
                      std::vector<Vector<Number>> &values) const;
 
  template <typename Number>
  void
  get_function_values(const ReadVector<Number> &fe_function,
                      const ArrayView<const types::global_dof_index> &indices,
                      std::vector<Number> &values) const;
 
  template <typename Number>
  void
  get_function_values(const ReadVector<Number> &fe_function,
                      const ArrayView<const types::global_dof_index> &indices,
                      std::vector<Vector<Number>> &values) const;
 
 
  template <typename Number>
  void
  get_function_values(const ReadVector<Number> &fe_function,
                      const ArrayView<const types::global_dof_index> &indices,
                      ArrayView<std::vector<Number>>                  values,
                      const bool quadrature_points_fastest) const;
 
 
  template <typename Number>
  void
  get_function_gradients(
    const ReadVector<Number>                 &fe_function,
    std::vector<Tensor<1, spacedim, Number>> &gradients) const;
 
  template <typename Number>
  void
  get_function_gradients(
    const ReadVector<Number>                              &fe_function,
    std::vector<std::vector<Tensor<1, spacedim, Number>>> &gradients) const;
 
  template <typename Number>
  void
  get_function_gradients(
    const ReadVector<Number>                       &fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<Tensor<1, spacedim, Number>>       &gradients) const;
 
  template <typename Number>
  void
  get_function_gradients(
    const ReadVector<Number>                           &fe_function,
    const ArrayView<const types::global_dof_index>     &indices,
    ArrayView<std::vector<Tensor<1, spacedim, Number>>> gradients,
    const bool quadrature_points_fastest = false) const;
 
 
  template <typename Number>
  void
  get_function_hessians(
    const ReadVector<Number>                 &fe_function,
    std::vector<Tensor<2, spacedim, Number>> &hessians) const;
 
  template <typename Number>
  void
  get_function_hessians(
    const ReadVector<Number>                              &fe_function,
    std::vector<std::vector<Tensor<2, spacedim, Number>>> &hessians,
    const bool quadrature_points_fastest = false) const;
 
  template <typename Number>
  void
  get_function_hessians(
    const ReadVector<Number>                       &fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<Tensor<2, spacedim, Number>>       &hessians) const;
 
  template <typename Number>
  void
  get_function_hessians(
    const ReadVector<Number>                           &fe_function,
    const ArrayView<const types::global_dof_index>     &indices,
    ArrayView<std::vector<Tensor<2, spacedim, Number>>> hessians,
    const bool quadrature_points_fastest = false) const;
 
  template <typename Number>
  void
  get_function_laplacians(const ReadVector<Number> &fe_function,
                          std::vector<Number>      &laplacians) const;
 
  template <typename Number>
  void
  get_function_laplacians(const ReadVector<Number>    &fe_function,
                          std::vector<Vector<Number>> &laplacians) const;
 
  template <typename Number>
  void
  get_function_laplacians(
    const ReadVector<Number>                       &fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<Number>                            &laplacians) const;
 
  template <typename Number>
  void
  get_function_laplacians(
    const ReadVector<Number>                       &fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<Vector<Number>>                    &laplacians) const;
 
  template <typename Number>
  void
  get_function_laplacians(
    const ReadVector<Number>                       &fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<std::vector<Number>>               &laplacians,
    const bool quadrature_points_fastest = false) const;
 
 
  template <typename Number>
  void
  get_function_third_derivatives(
    const ReadVector<Number>                 &fe_function,
    std::vector<Tensor<3, spacedim, Number>> &third_derivatives) const;
 
  template <typename Number>
  void
  get_function_third_derivatives(
    const ReadVector<Number>                              &fe_function,
    std::vector<std::vector<Tensor<3, spacedim, Number>>> &third_derivatives,
    const bool quadrature_points_fastest = false) const;
 
  template <typename Number>
  void
  get_function_third_derivatives(
    const ReadVector<Number>                       &fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<Tensor<3, spacedim, Number>>       &third_derivatives) const;
 
  template <typename Number>
  void
  get_function_third_derivatives(
    const ReadVector<Number>                           &fe_function,
    const ArrayView<const types::global_dof_index>     &indices,
    ArrayView<std::vector<Tensor<3, spacedim, Number>>> third_derivatives,
    const bool quadrature_points_fastest = false) const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  dof_indices() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  dof_indices_starting_at(const unsigned int start_dof_index) const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  dof_indices_ending_at(const unsigned int end_dof_index) const;
 
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  quadrature_point_indices() const;
 
  const Point<spacedim> &
  quadrature_point(const unsigned int q_point) const;
 
  const std::vector<Point<spacedim>> &
  get_quadrature_points() const;
 
  double
  JxW(const unsigned int q_point) const;
 
  const std::vector<double> &
  get_JxW_values() const;
 
  const DerivativeForm<1, dim, spacedim> &
  jacobian(const unsigned int q_point) const;
 
  const std::vector<DerivativeForm<1, dim, spacedim>> &
  get_jacobians() const;
 
  const DerivativeForm<2, dim, spacedim> &
  jacobian_grad(const unsigned int q_point) const;
 
  const std::vector<DerivativeForm<2, dim, spacedim>> &
  get_jacobian_grads() const;
 
  const Tensor<3, spacedim> &
  jacobian_pushed_forward_grad(const unsigned int q_point) const;
 
  const std::vector<Tensor<3, spacedim>> &
  get_jacobian_pushed_forward_grads() const;
 
  const DerivativeForm<3, dim, spacedim> &
  jacobian_2nd_derivative(const unsigned int q_point) const;
 
  const std::vector<DerivativeForm<3, dim, spacedim>> &
  get_jacobian_2nd_derivatives() const;
 
  const Tensor<4, spacedim> &
  jacobian_pushed_forward_2nd_derivative(const unsigned int q_point) const;
 
  const std::vector<Tensor<4, spacedim>> &
  get_jacobian_pushed_forward_2nd_derivatives() const;
 
  const DerivativeForm<4, dim, spacedim> &
  jacobian_3rd_derivative(const unsigned int q_point) const;
 
  const std::vector<DerivativeForm<4, dim, spacedim>> &
  get_jacobian_3rd_derivatives() const;
 
  const Tensor<5, spacedim> &
  jacobian_pushed_forward_3rd_derivative(const unsigned int q_point) const;
 
  const std::vector<Tensor<5, spacedim>> &
  get_jacobian_pushed_forward_3rd_derivatives() const;
 
  const DerivativeForm<1, spacedim, dim> &
  inverse_jacobian(const unsigned int q_point) const;
 
  const std::vector<DerivativeForm<1, spacedim, dim>> &
  get_inverse_jacobians() const;
 
  const Tensor<1, spacedim> &
  normal_vector(const unsigned int q_point) const;
 
  const std::vector<Tensor<1, spacedim>> &
  get_normal_vectors() const;
 
 
  const FEValuesViews::Scalar<dim, spacedim> &
  operator[](const FEValuesExtractors::Scalar &scalar) const;
 
  const FEValuesViews::Vector<dim, spacedim> &
  operator[](const FEValuesExtractors::Vector &vector) const;
 
  const FEValuesViews::SymmetricTensor<2, dim, spacedim> &
  operator[](const FEValuesExtractors::SymmetricTensor<2> &tensor) const;
 
 
  const FEValuesViews::Tensor<2, dim, spacedim> &
  operator[](const FEValuesExtractors::Tensor<2> &tensor) const;
 
 
  const Mapping<dim, spacedim> &
  get_mapping() const;
 
  const FiniteElement<dim, spacedim> &
  get_fe() const;
 
  UpdateFlags
  get_update_flags() const;
 
  typename Triangulation<dim, spacedim>::cell_iterator
  get_cell() const;
 
  CellSimilarity::Similarity
  get_cell_similarity() const;
 
  std::size_t
  memory_consumption() const;
  DeclException1(
    ExcAccessToUninitializedField,
    std::string,
    << "You are requesting information from an FEValues/FEFaceValues/FESubfaceValues "
    << "object for which this kind of information has not been computed. What "
    << "information these objects compute is determined by the update_* flags you "
    << "pass to the constructor. Here, the operation you are attempting requires "
    << "the <" << arg1
    << "> flag to be set, but it was apparently not specified "
    << "upon construction.");
 
  DeclExceptionMsg(ExcNotReinited,
                   "FEValues object is not reinit'ed to any cell");
 
  DeclExceptionMsg(
    ExcFEDontMatch,
    "The FiniteElement you provided to FEValues and the FiniteElement that belongs "
    "to the DoFHandler that provided the cell iterator do not match.");
  DeclException1(ExcShapeFunctionNotPrimitive,
                 int,
                 << "The shape function with index " << arg1
                 << " is not primitive, i.e. it is vector-valued and "
                 << "has more than one non-zero vector component. This "
                 << "function cannot be called for these shape functions. "
                 << "Maybe you want to use the same function with the "
                 << "_component suffix?");
 
  DeclExceptionMsg(
    ExcFENotPrimitive,
    "The given FiniteElement is not a primitive element but the requested operation "
    "only works for those. See FiniteElement::is_primitive() for more information.");
 
protected:
  class CellIteratorWrapper
  {
  public:
    DeclExceptionMsg(
      ExcNeedsDoFHandler,
      "You have previously called the FEValues::reinit() function with a "
      "cell iterator of type Triangulation<dim,spacedim>::cell_iterator. However, "
      "when you do this, you cannot call some functions in the FEValues "
      "class, such as the get_function_values/gradients/hessians/third_derivatives "
      "functions. If you need these functions, then you need to call "
      "FEValues::reinit() with an iterator type that allows to extract "
      "degrees of freedom, such as DoFHandler<dim,spacedim>::cell_iterator.");
 
    CellIteratorWrapper() = default;
 
    CellIteratorWrapper(
      const typename Triangulation<dim, spacedim>::cell_iterator &cell);
 
    CellIteratorWrapper(
      const typename DoFHandler<dim, spacedim>::cell_iterator &cell);
 
    CellIteratorWrapper(
      const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell);
 
    bool
    is_initialized() const;
 
    operator typename Triangulation<dim, spacedim>::cell_iterator() const;
 
    types::global_dof_index
    n_dofs_for_dof_handler() const;
 
    template <typename Number>
    void
    get_interpolated_dof_values(const ReadVector<Number> &in,
                                Vector<Number>           &out) const;
 
  private:
    std::optional<
      std::variant<typename Triangulation<dim, spacedim>::cell_iterator,
                   typename DoFHandler<dim, spacedim>::cell_iterator,
                   typename DoFHandler<dim, spacedim>::level_cell_iterator>>
      cell;
  };
 
  CellIteratorWrapper present_cell;
 
  boost::signals2::connection tria_listener_refinement;
 
  boost::signals2::connection tria_listener_mesh_transform;
 
  void
  invalidate_present_cell();
 
  void
  maybe_invalidate_previous_present_cell(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell);
 
  const SmartPointer<const Mapping<dim, spacedim>, FEValuesBase<dim, spacedim>>
    mapping;
 
  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
    mapping_data;
 
  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
    mapping_output;
 
  const SmartPointer<const FiniteElement<dim, spacedim>,
                     FEValuesBase<dim, spacedim>>
    fe;
 
  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
    fe_data;
 
  ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                     spacedim>
    finite_element_output;
 
 
  UpdateFlags update_flags;
 
  UpdateFlags
  compute_update_flags(const UpdateFlags update_flags) const;
 
  CellSimilarity::Similarity cell_similarity;
 
  void
  check_cell_similarity(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell);
 
private:
  ::internal::FEValuesViews::Cache<dim, spacedim> fe_values_views_cache;
 
  bool check_for_cell_similarity_allowed;
 
  // Make the view classes friends of this class, since they access internal
  // data.
  template <int, int>
  friend class FEValuesViews::Scalar;
  template <int, int>
  friend class FEValuesViews::Vector;
  template <int, int, int>
  friend class FEValuesViews::SymmetricTensor;
  template <int, int, int>
  friend class FEValuesViews::Tensor;
};
 
#ifndef DOXYGEN
 
/*---------------------- Inline functions: FEValuesBase ---------------------*/
 
template <int dim, int spacedim>
inline const FEValuesViews::Scalar<dim, spacedim> &
FEValuesBase<dim, spacedim>::operator[](
  const FEValuesExtractors::Scalar &scalar) const
{
  AssertIndexRange(scalar.component, fe_values_views_cache.scalars.size());
 
  return fe_values_views_cache.scalars[scalar.component].value_or_initialize(
    [scalar, this]() {
      return FEValuesViews::Scalar<dim, spacedim>(*this, scalar.component);
    });
}
 
 
 
template <int dim, int spacedim>
inline const FEValuesViews::Vector<dim, spacedim> &
FEValuesBase<dim, spacedim>::operator[](
  const FEValuesExtractors::Vector &vector) const
{
  AssertIndexRange(vector.first_vector_component,
                   fe_values_views_cache.vectors.size());
 
  return fe_values_views_cache.vectors[vector.first_vector_component]
    .value_or_initialize([vector, this]() {
      return FEValuesViews::Vector<dim, spacedim>(
        *this, vector.first_vector_component);
    });
}
 
 
 
template <int dim, int spacedim>
inline const FEValuesViews::SymmetricTensor<2, dim, spacedim> &
FEValuesBase<dim, spacedim>::operator[](
  const FEValuesExtractors::SymmetricTensor<2> &tensor) const
{
  Assert(
    tensor.first_tensor_component <
      fe_values_views_cache.symmetric_second_order_tensors.size(),
    ExcIndexRange(tensor.first_tensor_component,
                  0,
                  fe_values_views_cache.symmetric_second_order_tensors.size()));
 
  return fe_values_views_cache
    .symmetric_second_order_tensors[tensor.first_tensor_component]
    .value_or_initialize([tensor, this]() {
      return FEValuesViews::SymmetricTensor<2, dim, spacedim>(
        *this, tensor.first_tensor_component);
    });
}
 
 
 
template <int dim, int spacedim>
inline const FEValuesViews::Tensor<2, dim, spacedim> &
FEValuesBase<dim, spacedim>::operator[](
  const FEValuesExtractors::Tensor<2> &tensor) const
{
  AssertIndexRange(tensor.first_tensor_component,
                   fe_values_views_cache.second_order_tensors.size());
 
  return fe_values_views_cache
    .second_order_tensors[tensor.first_tensor_component]
    .value_or_initialize([tensor, this]() {
      return FEValuesViews::Tensor<2, dim, spacedim>(
        *this, tensor.first_tensor_component);
    });
}
 
 
 
template <int dim, int spacedim>
inline const double &
FEValuesBase<dim, spacedim>::shape_value(const unsigned int i,
                                         const unsigned int q_point) const
{
  AssertIndexRange(i, fe->n_dofs_per_cell());
  Assert(this->update_flags & update_values,
         ExcAccessToUninitializedField("update_values"));
  Assert(fe->is_primitive(i), ExcShapeFunctionNotPrimitive(i));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  // if the entire FE is primitive,
  // then we can take a short-cut:
  if (fe->is_primitive())
    return this->finite_element_output.shape_values(i, q_point);
  else
    {
      // otherwise, use the mapping
      // between shape function
      // numbers and rows. note that
      // by the assertions above, we
      // know that this particular
      // shape function is primitive,
      // so we can call
      // system_to_component_index
      const unsigned int row =
        this->finite_element_output
          .shape_function_to_row_table[i * fe->n_components() +
                                       fe->system_to_component_index(i).first];
      return this->finite_element_output.shape_values(row, q_point);
    }
}
 
 
 
template <int dim, int spacedim>
inline double
FEValuesBase<dim, spacedim>::shape_value_component(
  const unsigned int i,
  const unsigned int q_point,
  const unsigned int component) const
{
  AssertIndexRange(i, fe->n_dofs_per_cell());
  Assert(this->update_flags & update_values,
         ExcAccessToUninitializedField("update_values"));
  AssertIndexRange(component, fe->n_components());
  Assert(present_cell.is_initialized(), ExcNotReinited());
 
  // check whether the shape function
  // is non-zero at all within
  // this component:
  if (fe->get_nonzero_components(i)[component] == false)
    return 0;
 
  // look up the right row in the
  // table and take the data from
  // there
  const unsigned int row =
    this->finite_element_output
      .shape_function_to_row_table[i * fe->n_components() + component];
  return this->finite_element_output.shape_values(row, q_point);
}
 
 
 
template <int dim, int spacedim>
inline const Tensor<1, spacedim> &
FEValuesBase<dim, spacedim>::shape_grad(const unsigned int i,
                                        const unsigned int q_point) const
{
  AssertIndexRange(i, fe->n_dofs_per_cell());
  Assert(this->update_flags & update_gradients,
         ExcAccessToUninitializedField("update_gradients"));
  Assert(fe->is_primitive(i), ExcShapeFunctionNotPrimitive(i));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  // if the entire FE is primitive,
  // then we can take a short-cut:
  if (fe->is_primitive())
    return this->finite_element_output.shape_gradients[i][q_point];
  else
    {
      // otherwise, use the mapping
      // between shape function
      // numbers and rows. note that
      // by the assertions above, we
      // know that this particular
      // shape function is primitive,
      // so we can call
      // system_to_component_index
      const unsigned int row =
        this->finite_element_output
          .shape_function_to_row_table[i * fe->n_components() +
                                       fe->system_to_component_index(i).first];
      return this->finite_element_output.shape_gradients[row][q_point];
    }
}
 
 
 
template <int dim, int spacedim>
inline Tensor<1, spacedim>
FEValuesBase<dim, spacedim>::shape_grad_component(
  const unsigned int i,
  const unsigned int q_point,
  const unsigned int component) const
{
  AssertIndexRange(i, fe->n_dofs_per_cell());
  Assert(this->update_flags & update_gradients,
         ExcAccessToUninitializedField("update_gradients"));
  AssertIndexRange(component, fe->n_components());
  Assert(present_cell.is_initialized(), ExcNotReinited());
  // check whether the shape function
  // is non-zero at all within
  // this component:
  if (fe->get_nonzero_components(i)[component] == false)
    return Tensor<1, spacedim>();
 
  // look up the right row in the
  // table and take the data from
  // there
  const unsigned int row =
    this->finite_element_output
      .shape_function_to_row_table[i * fe->n_components() + component];
  return this->finite_element_output.shape_gradients[row][q_point];
}
 
 
 
template <int dim, int spacedim>
inline const Tensor<2, spacedim> &
FEValuesBase<dim, spacedim>::shape_hessian(const unsigned int i,
                                           const unsigned int q_point) const
{
  AssertIndexRange(i, fe->n_dofs_per_cell());
  Assert(this->update_flags & update_hessians,
         ExcAccessToUninitializedField("update_hessians"));
  Assert(fe->is_primitive(i), ExcShapeFunctionNotPrimitive(i));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  // if the entire FE is primitive,
  // then we can take a short-cut:
  if (fe->is_primitive())
    return this->finite_element_output.shape_hessians[i][q_point];
  else
    {
      // otherwise, use the mapping
      // between shape function
      // numbers and rows. note that
      // by the assertions above, we
      // know that this particular
      // shape function is primitive,
      // so we can call
      // system_to_component_index
      const unsigned int row =
        this->finite_element_output
          .shape_function_to_row_table[i * fe->n_components() +
                                       fe->system_to_component_index(i).first];
      return this->finite_element_output.shape_hessians[row][q_point];
    }
}
 
 
 
template <int dim, int spacedim>
inline Tensor<2, spacedim>
FEValuesBase<dim, spacedim>::shape_hessian_component(
  const unsigned int i,
  const unsigned int q_point,
  const unsigned int component) const
{
  AssertIndexRange(i, fe->n_dofs_per_cell());
  Assert(this->update_flags & update_hessians,
         ExcAccessToUninitializedField("update_hessians"));
  AssertIndexRange(component, fe->n_components());
  Assert(present_cell.is_initialized(), ExcNotReinited());
  // check whether the shape function
  // is non-zero at all within
  // this component:
  if (fe->get_nonzero_components(i)[component] == false)
    return Tensor<2, spacedim>();
 
  // look up the right row in the
  // table and take the data from
  // there
  const unsigned int row =
    this->finite_element_output
      .shape_function_to_row_table[i * fe->n_components() + component];
  return this->finite_element_output.shape_hessians[row][q_point];
}
 
 
 
template <int dim, int spacedim>
inline const Tensor<3, spacedim> &
FEValuesBase<dim, spacedim>::shape_3rd_derivative(
  const unsigned int i,
  const unsigned int q_point) const
{
  AssertIndexRange(i, fe->n_dofs_per_cell());
  Assert(this->update_flags & update_3rd_derivatives,
         ExcAccessToUninitializedField("update_3rd_derivatives"));
  Assert(fe->is_primitive(i), ExcShapeFunctionNotPrimitive(i));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  // if the entire FE is primitive,
  // then we can take a short-cut:
  if (fe->is_primitive())
    return this->finite_element_output.shape_3rd_derivatives[i][q_point];
  else
    {
      // otherwise, use the mapping
      // between shape function
      // numbers and rows. note that
      // by the assertions above, we
      // know that this particular
      // shape function is primitive,
      // so we can call
      // system_to_component_index
      const unsigned int row =
        this->finite_element_output
          .shape_function_to_row_table[i * fe->n_components() +
                                       fe->system_to_component_index(i).first];
      return this->finite_element_output.shape_3rd_derivatives[row][q_point];
    }
}
 
 
 
template <int dim, int spacedim>
inline Tensor<3, spacedim>
FEValuesBase<dim, spacedim>::shape_3rd_derivative_component(
  const unsigned int i,
  const unsigned int q_point,
  const unsigned int component) const
{
  AssertIndexRange(i, fe->n_dofs_per_cell());
  Assert(this->update_flags & update_3rd_derivatives,
         ExcAccessToUninitializedField("update_3rd_derivatives"));
  AssertIndexRange(component, fe->n_components());
  Assert(present_cell.is_initialized(), ExcNotReinited());
  // check whether the shape function
  // is non-zero at all within
  // this component:
  if (fe->get_nonzero_components(i)[component] == false)
    return Tensor<3, spacedim>();
 
  // look up the right row in the
  // table and take the data from
  // there
  const unsigned int row =
    this->finite_element_output
      .shape_function_to_row_table[i * fe->n_components() + component];
  return this->finite_element_output.shape_3rd_derivatives[row][q_point];
}
 
 
 
template <int dim, int spacedim>
inline const FiniteElement<dim, spacedim> &
FEValuesBase<dim, spacedim>::get_fe() const
{
  return *fe;
}
 
 
 
template <int dim, int spacedim>
inline const Mapping<dim, spacedim> &
FEValuesBase<dim, spacedim>::get_mapping() const
{
  return *mapping;
}
 
 
 
template <int dim, int spacedim>
inline UpdateFlags
FEValuesBase<dim, spacedim>::get_update_flags() const
{
  return this->update_flags;
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<Point<spacedim>> &
FEValuesBase<dim, spacedim>::get_quadrature_points() const
{
  Assert(this->update_flags & update_quadrature_points,
         ExcAccessToUninitializedField("update_quadrature_points"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.quadrature_points;
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<double> &
FEValuesBase<dim, spacedim>::get_JxW_values() const
{
  Assert(this->update_flags & update_JxW_values,
         ExcAccessToUninitializedField("update_JxW_values"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.JxW_values;
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<DerivativeForm<1, dim, spacedim>> &
FEValuesBase<dim, spacedim>::get_jacobians() const
{
  Assert(this->update_flags & update_jacobians,
         ExcAccessToUninitializedField("update_jacobians"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobians;
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<DerivativeForm<2, dim, spacedim>> &
FEValuesBase<dim, spacedim>::get_jacobian_grads() const
{
  Assert(this->update_flags & update_jacobian_grads,
         ExcAccessToUninitializedField("update_jacobians_grads"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_grads;
}
 
 
 
template <int dim, int spacedim>
inline const Tensor<3, spacedim> &
FEValuesBase<dim, spacedim>::jacobian_pushed_forward_grad(
  const unsigned int q_point) const
{
  Assert(this->update_flags & update_jacobian_pushed_forward_grads,
         ExcAccessToUninitializedField("update_jacobian_pushed_forward_grads"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_pushed_forward_grads[q_point];
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<Tensor<3, spacedim>> &
FEValuesBase<dim, spacedim>::get_jacobian_pushed_forward_grads() const
{
  Assert(this->update_flags & update_jacobian_pushed_forward_grads,
         ExcAccessToUninitializedField("update_jacobian_pushed_forward_grads"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_pushed_forward_grads;
}
 
 
 
template <int dim, int spacedim>
inline const DerivativeForm<3, dim, spacedim> &
FEValuesBase<dim, spacedim>::jacobian_2nd_derivative(
  const unsigned int q_point) const
{
  Assert(this->update_flags & update_jacobian_2nd_derivatives,
         ExcAccessToUninitializedField("update_jacobian_2nd_derivatives"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_2nd_derivatives[q_point];
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<DerivativeForm<3, dim, spacedim>> &
FEValuesBase<dim, spacedim>::get_jacobian_2nd_derivatives() const
{
  Assert(this->update_flags & update_jacobian_2nd_derivatives,
         ExcAccessToUninitializedField("update_jacobian_2nd_derivatives"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_2nd_derivatives;
}
 
 
 
template <int dim, int spacedim>
inline const Tensor<4, spacedim> &
FEValuesBase<dim, spacedim>::jacobian_pushed_forward_2nd_derivative(
  const unsigned int q_point) const
{
  Assert(this->update_flags & update_jacobian_pushed_forward_2nd_derivatives,
         ExcAccessToUninitializedField(
           "update_jacobian_pushed_forward_2nd_derivatives"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_pushed_forward_2nd_derivatives[q_point];
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<Tensor<4, spacedim>> &
FEValuesBase<dim, spacedim>::get_jacobian_pushed_forward_2nd_derivatives() const
{
  Assert(this->update_flags & update_jacobian_pushed_forward_2nd_derivatives,
         ExcAccessToUninitializedField(
           "update_jacobian_pushed_forward_2nd_derivatives"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_pushed_forward_2nd_derivatives;
}
 
 
 
template <int dim, int spacedim>
inline const DerivativeForm<4, dim, spacedim> &
FEValuesBase<dim, spacedim>::jacobian_3rd_derivative(
  const unsigned int q_point) const
{
  Assert(this->update_flags & update_jacobian_3rd_derivatives,
         ExcAccessToUninitializedField("update_jacobian_3rd_derivatives"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_3rd_derivatives[q_point];
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<DerivativeForm<4, dim, spacedim>> &
FEValuesBase<dim, spacedim>::get_jacobian_3rd_derivatives() const
{
  Assert(this->update_flags & update_jacobian_3rd_derivatives,
         ExcAccessToUninitializedField("update_jacobian_3rd_derivatives"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_3rd_derivatives;
}
 
 
 
template <int dim, int spacedim>
inline const Tensor<5, spacedim> &
FEValuesBase<dim, spacedim>::jacobian_pushed_forward_3rd_derivative(
  const unsigned int q_point) const
{
  Assert(this->update_flags & update_jacobian_pushed_forward_3rd_derivatives,
         ExcAccessToUninitializedField(
           "update_jacobian_pushed_forward_3rd_derivatives"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_pushed_forward_3rd_derivatives[q_point];
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<Tensor<5, spacedim>> &
FEValuesBase<dim, spacedim>::get_jacobian_pushed_forward_3rd_derivatives() const
{
  Assert(this->update_flags & update_jacobian_pushed_forward_3rd_derivatives,
         ExcAccessToUninitializedField(
           "update_jacobian_pushed_forward_3rd_derivatives"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.jacobian_pushed_forward_3rd_derivatives;
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<DerivativeForm<1, spacedim, dim>> &
FEValuesBase<dim, spacedim>::get_inverse_jacobians() const
{
  Assert(this->update_flags & update_inverse_jacobians,
         ExcAccessToUninitializedField("update_inverse_jacobians"));
  Assert(present_cell.is_initialized(), ExcNotReinited());
  return this->mapping_output.inverse_jacobians;
}
 
 
 
template <int dim, int spacedim>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FEValuesBase<dim, spacedim>::dof_indices() const
{
  return {0U, dofs_per_cell};
}
 
 
 
template <int dim, int spacedim>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FEValuesBase<dim, spacedim>::dof_indices_starting_at(
  const unsigned int start_dof_index) const
{
  Assert(start_dof_index <= dofs_per_cell,
         ExcIndexRange(start_dof_index, 0, dofs_per_cell + 1));
  return {start_dof_index, dofs_per_cell};
}
 
 
 
template <int dim, int spacedim>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FEValuesBase<dim, spacedim>::dof_indices_ending_at(
  const unsigned int end_dof_index) const
{
  Assert(end_dof_index < dofs_per_cell,
         ExcIndexRange(end_dof_index, 0, dofs_per_cell));
  return {0U, end_dof_index + 1};
}
 
 
 
template <int dim, int spacedim>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FEValuesBase<dim, spacedim>::quadrature_point_indices() const
{
  return {0U, n_quadrature_points};
}
 
 
 
template <int dim, int spacedim>
inline const Point<spacedim> &
FEValuesBase<dim, spacedim>::quadrature_point(const unsigned int q_point) const
{
  Assert(this->update_flags & update_quadrature_points,
         ExcAccessToUninitializedField("update_quadrature_points"));
  AssertIndexRange(q_point, this->mapping_output.quadrature_points.size());
  Assert(present_cell.is_initialized(), ExcNotReinited());
 
  return this->mapping_output.quadrature_points[q_point];
}
 
 
 
template <int dim, int spacedim>
inline double
FEValuesBase<dim, spacedim>::JxW(const unsigned int q_point) const
{
  Assert(this->update_flags & update_JxW_values,
         ExcAccessToUninitializedField("update_JxW_values"));
  AssertIndexRange(q_point, this->mapping_output.JxW_values.size());
  Assert(present_cell.is_initialized(), ExcNotReinited());
 
  return this->mapping_output.JxW_values[q_point];
}
 
 
 
template <int dim, int spacedim>
inline const DerivativeForm<1, dim, spacedim> &
FEValuesBase<dim, spacedim>::jacobian(const unsigned int q_point) const
{
  Assert(this->update_flags & update_jacobians,
         ExcAccessToUninitializedField("update_jacobians"));
  AssertIndexRange(q_point, this->mapping_output.jacobians.size());
  Assert(present_cell.is_initialized(), ExcNotReinited());
 
  return this->mapping_output.jacobians[q_point];
}
 
 
 
template <int dim, int spacedim>
inline const DerivativeForm<2, dim, spacedim> &
FEValuesBase<dim, spacedim>::jacobian_grad(const unsigned int q_point) const
{
  Assert(this->update_flags & update_jacobian_grads,
         ExcAccessToUninitializedField("update_jacobians_grads"));
  AssertIndexRange(q_point, this->mapping_output.jacobian_grads.size());
  Assert(present_cell.is_initialized(), ExcNotReinited());
 
  return this->mapping_output.jacobian_grads[q_point];
}
 
 
 
template <int dim, int spacedim>
inline const DerivativeForm<1, spacedim, dim> &
FEValuesBase<dim, spacedim>::inverse_jacobian(const unsigned int q_point) const
{
  Assert(this->update_flags & update_inverse_jacobians,
         ExcAccessToUninitializedField("update_inverse_jacobians"));
  AssertIndexRange(q_point, this->mapping_output.inverse_jacobians.size());
  Assert(present_cell.is_initialized(), ExcNotReinited());
 
  return this->mapping_output.inverse_jacobians[q_point];
}
 
 
 
template <int dim, int spacedim>
inline const Tensor<1, spacedim> &
FEValuesBase<dim, spacedim>::normal_vector(const unsigned int q_point) const
{
  Assert(this->update_flags & update_normal_vectors,
         (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
           "update_normal_vectors")));
  AssertIndexRange(q_point, this->mapping_output.normal_vectors.size());
  Assert(present_cell.is_initialized(), ExcNotReinited());
 
  return this->mapping_output.normal_vectors[q_point];
}
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif