fe_values.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, 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.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_fe_values_h
#define dealii_fe_values_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/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 <algorithm>
#include <memory>
#include <type_traits>
 
 
// dummy include in order to have the
// definition of PetscScalar available
// without including other PETSc stuff
#ifdef DEAL_II_WITH_PETSC
#  include <petsc.h>
#endif
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declaration
#ifndef DOXYGEN
template <int dim, int spacedim = dim>
class FEValuesBase;
#endif
 
namespace internal
{
  template <int dim, class NumberType = double>
  struct CurlType;
 
  template <class NumberType>
  struct CurlType<1, NumberType>
  {
    using type = Tensor<1, 1, NumberType>;
  };
 
  template <class NumberType>
  struct CurlType<2, NumberType>
  {
    using type = Tensor<1, 1, NumberType>;
  };
 
  template <class NumberType>
  struct CurlType<3, NumberType>
  {
    using type = Tensor<1, 3, NumberType>;
  };
} // namespace internal
 
 
 
namespace FEValuesViews
{
  template <int dim, int spacedim = dim>
  class Scalar
  {
  public:
    using value_type = double;
 
    using gradient_type = ::Tensor<1, spacedim>;
 
    using hessian_type = ::Tensor<2, spacedim>;
 
    using third_derivative_type = ::Tensor<3, spacedim>;
 
    template <typename Number>
    using solution_value_type = typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_gradient_type =
      typename ProductType<Number, gradient_type>::type;
 
    template <typename Number>
    using solution_laplacian_type =
      typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_hessian_type =
      typename ProductType<Number, hessian_type>::type;
 
    template <typename Number>
    using solution_third_derivative_type =
      typename ProductType<Number, third_derivative_type>::type;
 
    template <typename Number>
    struct DEAL_II_DEPRECATED OutputType
    {
      using value_type =
        typename ProductType<Number,
                             typename Scalar<dim, spacedim>::value_type>::type;
 
      using gradient_type = typename ProductType<
        Number,
        typename Scalar<dim, spacedim>::gradient_type>::type;
 
      using laplacian_type =
        typename ProductType<Number,
                             typename Scalar<dim, spacedim>::value_type>::type;
 
      using hessian_type = typename ProductType<
        Number,
        typename Scalar<dim, spacedim>::hessian_type>::type;
 
      using third_derivative_type = typename ProductType<
        Number,
        typename Scalar<dim, spacedim>::third_derivative_type>::type;
    };
 
    struct ShapeFunctionData
    {
      bool is_nonzero_shape_function_component;
 
      unsigned int row_index;
    };
 
    Scalar();
 
    Scalar(const FEValuesBase<dim, spacedim> &fe_values_base,
           const unsigned int                 component);
 
    Scalar(const Scalar<dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    Scalar(Scalar<dim, spacedim> &&) = default;
 
    ~Scalar() = default;
 
    Scalar &
    operator=(const Scalar<dim, spacedim> &) = delete;
 
    Scalar &
    operator=(Scalar<dim, spacedim> &&) noexcept = default;
 
    value_type
    value(const unsigned int shape_function, const unsigned int q_point) const;
 
    gradient_type
    gradient(const unsigned int shape_function,
             const unsigned int q_point) const;
 
    hessian_type
    hessian(const unsigned int shape_function,
            const unsigned int q_point) const;
 
    third_derivative_type
    third_derivative(const unsigned int shape_function,
                     const unsigned int q_point) const;
 
    template <class InputVector>
    void
    get_function_values(
      const InputVector &fe_function,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <class InputVector>
    void
    get_function_values_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <class InputVector>
    void
    get_function_gradients(
      const InputVector &fe_function,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
    template <class InputVector>
    void
    get_function_gradients_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
    template <class InputVector>
    void
    get_function_hessians(
      const InputVector &fe_function,
      std::vector<solution_hessian_type<typename InputVector::value_type>>
        &hessians) const;
 
    template <class InputVector>
    void
    get_function_hessians_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_hessian_type<typename InputVector::value_type>>
        &hessians) const;
 
 
    template <class InputVector>
    void
    get_function_laplacians(
      const InputVector &fe_function,
      std::vector<solution_laplacian_type<typename InputVector::value_type>>
        &laplacians) const;
 
    template <class InputVector>
    void
    get_function_laplacians_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_laplacian_type<typename InputVector::value_type>>
        &laplacians) const;
 
 
    template <class InputVector>
    void
    get_function_third_derivatives(
      const InputVector &fe_function,
      std::vector<
        solution_third_derivative_type<typename InputVector::value_type>>
        &third_derivatives) const;
 
    template <class InputVector>
    void
    get_function_third_derivatives_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<
        solution_third_derivative_type<typename InputVector::value_type>>
        &third_derivatives) const;
 
 
  private:
    const SmartPointer<const FEValuesBase<dim, spacedim>> fe_values;
 
    const unsigned int component;
 
    std::vector<ShapeFunctionData> shape_function_data;
  };
 
 
 
  template <int dim, int spacedim = dim>
  class Vector
  {
  public:
    using value_type = ::Tensor<1, spacedim>;
 
    using gradient_type = ::Tensor<2, spacedim>;
 
    using symmetric_gradient_type = ::SymmetricTensor<2, spacedim>;
 
    using divergence_type = double;
 
    using curl_type = typename ::internal::CurlType<spacedim>::type;
 
    using hessian_type = ::Tensor<3, spacedim>;
 
    using third_derivative_type = ::Tensor<4, spacedim>;
 
    template <typename Number>
    using solution_value_type = typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_gradient_type =
      typename ProductType<Number, gradient_type>::type;
 
    template <typename Number>
    using solution_symmetric_gradient_type =
      typename ProductType<Number, symmetric_gradient_type>::type;
 
    template <typename Number>
    using solution_divergence_type =
      typename ProductType<Number, divergence_type>::type;
 
    template <typename Number>
    using solution_laplacian_type =
      typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_curl_type = typename ProductType<Number, curl_type>::type;
 
    template <typename Number>
    using solution_hessian_type =
      typename ProductType<Number, hessian_type>::type;
 
    template <typename Number>
    using solution_third_derivative_type =
      typename ProductType<Number, third_derivative_type>::type;
 
    template <typename Number>
    struct DEAL_II_DEPRECATED OutputType
    {
      using value_type =
        typename ProductType<Number,
                             typename Vector<dim, spacedim>::value_type>::type;
 
      using gradient_type = typename ProductType<
        Number,
        typename Vector<dim, spacedim>::gradient_type>::type;
 
      using symmetric_gradient_type = typename ProductType<
        Number,
        typename Vector<dim, spacedim>::symmetric_gradient_type>::type;
 
      using divergence_type = typename ProductType<
        Number,
        typename Vector<dim, spacedim>::divergence_type>::type;
 
      using laplacian_type =
        typename ProductType<Number,
                             typename Vector<dim, spacedim>::value_type>::type;
 
      using curl_type =
        typename ProductType<Number,
                             typename Vector<dim, spacedim>::curl_type>::type;
 
      using hessian_type = typename ProductType<
        Number,
        typename Vector<dim, spacedim>::hessian_type>::type;
 
      using third_derivative_type = typename ProductType<
        Number,
        typename Vector<dim, spacedim>::third_derivative_type>::type;
    };
 
    struct ShapeFunctionData
    {
      bool is_nonzero_shape_function_component[spacedim];
 
      unsigned int row_index[spacedim];
 
      int          single_nonzero_component;
      unsigned int single_nonzero_component_index;
    };
 
    Vector();
 
    Vector(const FEValuesBase<dim, spacedim> &fe_values_base,
           const unsigned int                 first_vector_component);
 
    Vector(const Vector<dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    Vector(Vector<dim, spacedim> &&) = default;
 
    ~Vector() = default;
 
    Vector &
    operator=(const Vector<dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    Vector &
    operator=(Vector<dim, spacedim> &&) = default; // NOLINT
 
    value_type
    value(const unsigned int shape_function, const unsigned int q_point) const;
 
    gradient_type
    gradient(const unsigned int shape_function,
             const unsigned int q_point) const;
 
    symmetric_gradient_type
    symmetric_gradient(const unsigned int shape_function,
                       const unsigned int q_point) const;
 
    divergence_type
    divergence(const unsigned int shape_function,
               const unsigned int q_point) const;
 
    curl_type
    curl(const unsigned int shape_function, const unsigned int q_point) const;
 
    hessian_type
    hessian(const unsigned int shape_function,
            const unsigned int q_point) const;
 
    third_derivative_type
    third_derivative(const unsigned int shape_function,
                     const unsigned int q_point) const;
 
    template <class InputVector>
    void
    get_function_values(
      const InputVector &fe_function,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <class InputVector>
    void
    get_function_values_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <class InputVector>
    void
    get_function_gradients(
      const InputVector &fe_function,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
    template <class InputVector>
    void
    get_function_gradients_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
    template <class InputVector>
    void
    get_function_symmetric_gradients(
      const InputVector &fe_function,
      std::vector<
        solution_symmetric_gradient_type<typename InputVector::value_type>>
        &symmetric_gradients) const;
 
    template <class InputVector>
    void
    get_function_symmetric_gradients_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<
        solution_symmetric_gradient_type<typename InputVector::value_type>>
        &symmetric_gradients) const;
 
    template <class InputVector>
    void
    get_function_divergences(
      const InputVector &fe_function,
      std::vector<solution_divergence_type<typename InputVector::value_type>>
        &divergences) const;
 
    template <class InputVector>
    void
    get_function_divergences_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_divergence_type<typename InputVector::value_type>>
        &divergences) const;
 
    template <class InputVector>
    void
    get_function_curls(
      const InputVector &fe_function,
      std::vector<solution_curl_type<typename InputVector::value_type>> &curls)
      const;
 
    template <class InputVector>
    void
    get_function_curls_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_curl_type<typename InputVector::value_type>> &curls)
      const;
 
    template <class InputVector>
    void
    get_function_hessians(
      const InputVector &fe_function,
      std::vector<solution_hessian_type<typename InputVector::value_type>>
        &hessians) const;
 
    template <class InputVector>
    void
    get_function_hessians_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_hessian_type<typename InputVector::value_type>>
        &hessians) const;
 
    template <class InputVector>
    void
    get_function_laplacians(
      const InputVector &fe_function,
      std::vector<solution_laplacian_type<typename InputVector::value_type>>
        &laplacians) const;
 
    template <class InputVector>
    void
    get_function_laplacians_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_laplacian_type<typename InputVector::value_type>>
        &laplacians) const;
 
    template <class InputVector>
    void
    get_function_third_derivatives(
      const InputVector &fe_function,
      std::vector<
        solution_third_derivative_type<typename InputVector::value_type>>
        &third_derivatives) const;
 
    template <class InputVector>
    void
    get_function_third_derivatives_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<
        solution_third_derivative_type<typename InputVector::value_type>>
        &third_derivatives) const;
 
  private:
    const SmartPointer<const FEValuesBase<dim, spacedim>> fe_values;
 
    const unsigned int first_vector_component;
 
    std::vector<ShapeFunctionData> shape_function_data;
  };
 
 
  template <int rank, int dim, int spacedim = dim>
  class SymmetricTensor;
 
  template <int dim, int spacedim>
  class SymmetricTensor<2, dim, spacedim>
  {
  public:
    using value_type = ::SymmetricTensor<2, spacedim>;
 
    using divergence_type = ::Tensor<1, spacedim>;
 
    template <typename Number>
    using solution_value_type = typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_divergence_type =
      typename ProductType<Number, divergence_type>::type;
 
 
    template <typename Number>
    struct DEAL_II_DEPRECATED OutputType
    {
      using value_type = typename ProductType<
        Number,
        typename SymmetricTensor<2, dim, spacedim>::value_type>::type;
 
      using divergence_type = typename ProductType<
        Number,
        typename SymmetricTensor<2, dim, spacedim>::divergence_type>::type;
    };
 
    struct ShapeFunctionData
    {
      bool is_nonzero_shape_function_component
        [value_type::n_independent_components];
 
      unsigned int row_index[value_type::n_independent_components];
 
      int single_nonzero_component;
 
      unsigned int single_nonzero_component_index;
    };
 
    SymmetricTensor();
 
    SymmetricTensor(const FEValuesBase<dim, spacedim> &fe_values_base,
                    const unsigned int                 first_tensor_component);
 
    SymmetricTensor(const SymmetricTensor<2, dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    SymmetricTensor(SymmetricTensor<2, dim, spacedim> &&) = default;
 
    SymmetricTensor &
    operator=(const SymmetricTensor<2, dim, spacedim> &) = delete;
 
    SymmetricTensor &
    operator=(SymmetricTensor<2, dim, spacedim> &&) noexcept = default;
 
    value_type
    value(const unsigned int shape_function, const unsigned int q_point) const;
 
    divergence_type
    divergence(const unsigned int shape_function,
               const unsigned int q_point) const;
 
    template <class InputVector>
    void
    get_function_values(
      const InputVector &fe_function,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <class InputVector>
    void
    get_function_values_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <class InputVector>
    void
    get_function_divergences(
      const InputVector &fe_function,
      std::vector<solution_divergence_type<typename InputVector::value_type>>
        &divergences) const;
 
    template <class InputVector>
    void
    get_function_divergences_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_divergence_type<typename InputVector::value_type>>
        &divergences) const;
 
  private:
    const SmartPointer<const FEValuesBase<dim, spacedim>> fe_values;
 
    const unsigned int first_tensor_component;
 
    std::vector<ShapeFunctionData> shape_function_data;
  };
 
 
  template <int rank, int dim, int spacedim = dim>
  class Tensor;
 
  template <int dim, int spacedim>
  class Tensor<2, dim, spacedim>
  {
  public:
    using value_type = ::Tensor<2, spacedim>;
 
    using divergence_type = ::Tensor<1, spacedim>;
 
    using gradient_type = ::Tensor<3, spacedim>;
 
    template <typename Number>
    using solution_value_type = typename ProductType<Number, value_type>::type;
 
    template <typename Number>
    using solution_divergence_type =
      typename ProductType<Number, divergence_type>::type;
 
    template <typename Number>
    using solution_gradient_type =
      typename ProductType<Number, gradient_type>::type;
 
 
    template <typename Number>
    struct DEAL_II_DEPRECATED OutputType
    {
      using value_type = typename ProductType<
        Number,
        typename Tensor<2, dim, spacedim>::value_type>::type;
 
      using divergence_type = typename ProductType<
        Number,
        typename Tensor<2, dim, spacedim>::divergence_type>::type;
 
      using gradient_type = typename ProductType<
        Number,
        typename Tensor<2, dim, spacedim>::gradient_type>::type;
    };
 
    struct ShapeFunctionData
    {
      bool is_nonzero_shape_function_component
        [value_type::n_independent_components];
 
      unsigned int row_index[value_type::n_independent_components];
 
      int single_nonzero_component;
 
      unsigned int single_nonzero_component_index;
    };
 
    Tensor();
 
    Tensor(const Tensor<2, dim, spacedim> &) = delete;
 
    // NOLINTNEXTLINE OSX does not compile with noexcept
    Tensor(Tensor<2, dim, spacedim> &&) = default;
 
    ~Tensor() = default;
 
    Tensor(const FEValuesBase<dim, spacedim> &fe_values_base,
           const unsigned int                 first_tensor_component);
 
 
    Tensor &
    operator=(const Tensor<2, dim, spacedim> &) = delete;
 
    Tensor &
    operator=(Tensor<2, dim, spacedim> &&) = default; // NOLINT
 
    value_type
    value(const unsigned int shape_function, const unsigned int q_point) const;
 
    divergence_type
    divergence(const unsigned int shape_function,
               const unsigned int q_point) const;
 
    gradient_type
    gradient(const unsigned int shape_function,
             const unsigned int q_point) const;
 
    template <class InputVector>
    void
    get_function_values(
      const InputVector &fe_function,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <class InputVector>
    void
    get_function_values_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <class InputVector>
    void
    get_function_divergences(
      const InputVector &fe_function,
      std::vector<solution_divergence_type<typename InputVector::value_type>>
        &divergences) const;
 
    template <class InputVector>
    void
    get_function_divergences_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_divergence_type<typename InputVector::value_type>>
        &divergences) const;
 
    template <class InputVector>
    void
    get_function_gradients(
      const InputVector &fe_function,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
    template <class InputVector>
    void
    get_function_gradients_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
  private:
    const SmartPointer<const FEValuesBase<dim, spacedim>> fe_values;
 
    const unsigned int first_tensor_component;
 
    std::vector<ShapeFunctionData> shape_function_data;
  };
 
} // namespace FEValuesViews
 
 
namespace internal
{
  namespace FEValuesViews
  {
    template <int dim, int spacedim, typename Extractor>
    struct ViewType
    {};
 
    template <int dim, int spacedim>
    struct ViewType<dim, spacedim, FEValuesExtractors::Scalar>
    {
      using type = typename ::FEValuesViews::Scalar<dim, spacedim>;
    };
 
    template <int dim, int spacedim>
    struct ViewType<dim, spacedim, FEValuesExtractors::Vector>
    {
      using type = typename ::FEValuesViews::Vector<dim, spacedim>;
    };
 
    template <int dim, int spacedim, int rank>
    struct ViewType<dim, spacedim, FEValuesExtractors::Tensor<rank>>
    {
      using type = typename ::FEValuesViews::Tensor<rank, dim, spacedim>;
    };
 
    template <int dim, int spacedim, int rank>
    struct ViewType<dim, spacedim, FEValuesExtractors::SymmetricTensor<rank>>
    {
      using type =
        typename ::FEValuesViews::SymmetricTensor<rank, dim, spacedim>;
    };
 
    template <int dim, int spacedim>
    struct Cache
    {
      std::vector<::FEValuesViews::Scalar<dim, spacedim>> scalars;
      std::vector<::FEValuesViews::Vector<dim, spacedim>> vectors;
      std::vector<::FEValuesViews::SymmetricTensor<2, dim, spacedim>>
        symmetric_second_order_tensors;
      std::vector<::FEValuesViews::Tensor<2, dim, spacedim>>
        second_order_tensors;
 
      Cache(const FEValuesBase<dim, spacedim> &fe_values);
    };
  } // namespace FEValuesViews
} // namespace internal
 
namespace FEValuesViews
{
  template <int dim, int spacedim, typename Extractor>
  using View = typename ::internal::FEValuesViews::
    ViewType<dim, spacedim, Extractor>::type;
} // namespace FEValuesViews
 
 
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;
 
 
 
  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 <class InputVector>
  void
  get_function_values(
    const InputVector &                            fe_function,
    std::vector<typename InputVector::value_type> &values) const;
 
  template <class InputVector>
  void
  get_function_values(
    const InputVector &                                    fe_function,
    std::vector<Vector<typename InputVector::value_type>> &values) const;
 
  template <class InputVector>
  void
  get_function_values(
    const InputVector &                             fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<typename InputVector::value_type> & values) const;
 
  template <class InputVector>
  void
  get_function_values(
    const InputVector &                                    fe_function,
    const ArrayView<const types::global_dof_index> &       indices,
    std::vector<Vector<typename InputVector::value_type>> &values) const;
 
 
  template <class InputVector>
  void
  get_function_values(
    const InputVector &                                      fe_function,
    const ArrayView<const types::global_dof_index> &         indices,
    ArrayView<std::vector<typename InputVector::value_type>> values,
    const bool quadrature_points_fastest) const;
 
 
  template <class InputVector>
  void
  get_function_gradients(
    const InputVector &fe_function,
    std::vector<Tensor<1, spacedim, typename InputVector::value_type>>
      &gradients) const;
 
  template <class InputVector>
  void
  get_function_gradients(
    const InputVector &fe_function,
    std::vector<
      std::vector<Tensor<1, spacedim, typename InputVector::value_type>>>
      &gradients) const;
 
  template <class InputVector>
  void
  get_function_gradients(
    const InputVector &                             fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<Tensor<1, spacedim, typename InputVector::value_type>>
      &gradients) const;
 
  template <class InputVector>
  void
  get_function_gradients(
    const InputVector &                             fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    ArrayView<
      std::vector<Tensor<1, spacedim, typename InputVector::value_type>>>
               gradients,
    const bool quadrature_points_fastest = false) const;
 
 
  template <class InputVector>
  void
  get_function_hessians(
    const InputVector &fe_function,
    std::vector<Tensor<2, spacedim, typename InputVector::value_type>>
      &hessians) const;
 
  template <class InputVector>
  void
  get_function_hessians(
    const InputVector &fe_function,
    std::vector<
      std::vector<Tensor<2, spacedim, typename InputVector::value_type>>>
      &        hessians,
    const bool quadrature_points_fastest = false) const;
 
  template <class InputVector>
  void
  get_function_hessians(
    const InputVector &                             fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<Tensor<2, spacedim, typename InputVector::value_type>>
      &hessians) const;
 
  template <class InputVector>
  void
  get_function_hessians(
    const InputVector &                             fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    ArrayView<
      std::vector<Tensor<2, spacedim, typename InputVector::value_type>>>
               hessians,
    const bool quadrature_points_fastest = false) const;
 
  template <class InputVector>
  void
  get_function_laplacians(
    const InputVector &                            fe_function,
    std::vector<typename InputVector::value_type> &laplacians) const;
 
  template <class InputVector>
  void
  get_function_laplacians(
    const InputVector &                                    fe_function,
    std::vector<Vector<typename InputVector::value_type>> &laplacians) const;
 
  template <class InputVector>
  void
  get_function_laplacians(
    const InputVector &                             fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<typename InputVector::value_type> & laplacians) const;
 
  template <class InputVector>
  void
  get_function_laplacians(
    const InputVector &                                    fe_function,
    const ArrayView<const types::global_dof_index> &       indices,
    std::vector<Vector<typename InputVector::value_type>> &laplacians) const;
 
  template <class InputVector>
  void
  get_function_laplacians(
    const InputVector &                                         fe_function,
    const ArrayView<const types::global_dof_index> &            indices,
    std::vector<std::vector<typename InputVector::value_type>> &laplacians,
    const bool quadrature_points_fastest = false) const;
 
 
  template <class InputVector>
  void
  get_function_third_derivatives(
    const InputVector &fe_function,
    std::vector<Tensor<3, spacedim, typename InputVector::value_type>>
      &third_derivatives) const;
 
  template <class InputVector>
  void
  get_function_third_derivatives(
    const InputVector &fe_function,
    std::vector<
      std::vector<Tensor<3, spacedim, typename InputVector::value_type>>>
      &        third_derivatives,
    const bool quadrature_points_fastest = false) const;
 
  template <class InputVector>
  void
  get_function_third_derivatives(
    const InputVector &                             fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    std::vector<Tensor<3, spacedim, typename InputVector::value_type>>
      &third_derivatives) const;
 
  template <class InputVector>
  void
  get_function_third_derivatives(
    const InputVector &                             fe_function,
    const ArrayView<const types::global_dof_index> &indices,
    ArrayView<
      std::vector<Tensor<3, spacedim, typename InputVector::value_type>>>
               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;
 
  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 CellIteratorContainer
  {
  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.");
 
    CellIteratorContainer();
 
    template <bool lda>
    CellIteratorContainer(
      const TriaIterator<DoFCellAccessor<dim, spacedim, lda>> &cell);
 
    CellIteratorContainer(
      const typename Triangulation<dim, spacedim>::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 VectorType>
    void
    get_interpolated_dof_values(
      const VectorType &                       in,
      Vector<typename VectorType::value_type> &out) const;
 
    void
    get_interpolated_dof_values(const IndexSet &              in,
                                Vector<IndexSet::value_type> &out) const;
 
  private:
    bool                                                 initialized;
    typename Triangulation<dim, spacedim>::cell_iterator cell;
    const DoFHandler<dim, spacedim> *                    dof_handler;
    bool                                                 level_dof_access;
  };
 
  CellIteratorContainer 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;
 
  // 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;
};
 
 
 
template <int dim, int spacedim = dim>
class FEValues : public FEValuesBase<dim, spacedim>
{
public:
  static constexpr unsigned int integral_dimension = dim;
 
  FEValues(const Mapping<dim, spacedim> &      mapping,
           const FiniteElement<dim, spacedim> &fe,
           const Quadrature<dim> &             quadrature,
           const UpdateFlags                   update_flags);
 
  FEValues(const Mapping<dim, spacedim> &      mapping,
           const FiniteElement<dim, spacedim> &fe,
           const hp::QCollection<dim> &        quadrature,
           const UpdateFlags                   update_flags);
 
  FEValues(const FiniteElement<dim, spacedim> &fe,
           const Quadrature<dim> &             quadrature,
           const UpdateFlags                   update_flags);
 
  FEValues(const FiniteElement<dim, spacedim> &fe,
           const hp::QCollection<dim> &        quadrature,
           const UpdateFlags                   update_flags);
 
  template <bool level_dof_access>
  void
  reinit(
    const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell);
 
  void
  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell);
 
  const Quadrature<dim> &
  get_quadrature() const;
 
  std::size_t
  memory_consumption() const;
 
  const FEValues<dim, spacedim> &
  get_present_fe_values() const;
 
private:
  const Quadrature<dim> quadrature;
 
  void
  initialize(const UpdateFlags update_flags);
 
  void
  do_reinit();
};
 
 
template <int dim, int spacedim = dim>
class FEFaceValuesBase : public FEValuesBase<dim, spacedim>
{
public:
  static constexpr unsigned int integral_dimension = dim - 1;
 
  FEFaceValuesBase(const unsigned int                  dofs_per_cell,
                   const UpdateFlags                   update_flags,
                   const Mapping<dim, spacedim> &      mapping,
                   const FiniteElement<dim, spacedim> &fe,
                   const Quadrature<dim - 1> &         quadrature);
 
  FEFaceValuesBase(const unsigned int                  dofs_per_cell,
                   const UpdateFlags                   update_flags,
                   const Mapping<dim, spacedim> &      mapping,
                   const FiniteElement<dim, spacedim> &fe,
                   const hp::QCollection<dim - 1> &    quadrature);
 
  const Tensor<1, spacedim> &
  boundary_form(const unsigned int q_point) const;
 
  const std::vector<Tensor<1, spacedim>> &
  get_boundary_forms() const;
 
  unsigned int
  get_face_number() const;
 
  unsigned int
  get_face_index() const;
 
  const Quadrature<dim - 1> &
  get_quadrature() const;
 
  std::size_t
  memory_consumption() const;
 
protected:
  unsigned int present_face_no;
 
  unsigned int present_face_index;
 
  const hp::QCollection<dim - 1> quadrature;
};
 
 
 
template <int dim, int spacedim = dim>
class FEFaceValues : public FEFaceValuesBase<dim, spacedim>
{
public:
  static constexpr unsigned int dimension = dim;
 
  static constexpr unsigned int space_dimension = spacedim;
 
  static constexpr unsigned int integral_dimension = dim - 1;
 
  FEFaceValues(const Mapping<dim, spacedim> &      mapping,
               const FiniteElement<dim, spacedim> &fe,
               const Quadrature<dim - 1> &         quadrature,
               const UpdateFlags                   update_flags);
 
  FEFaceValues(const Mapping<dim, spacedim> &      mapping,
               const FiniteElement<dim, spacedim> &fe,
               const hp::QCollection<dim - 1> &    quadrature,
               const UpdateFlags                   update_flags);
 
  FEFaceValues(const FiniteElement<dim, spacedim> &fe,
               const Quadrature<dim - 1> &         quadrature,
               const UpdateFlags                   update_flags);
 
  FEFaceValues(const FiniteElement<dim, spacedim> &fe,
               const hp::QCollection<dim - 1> &    quadrature,
               const UpdateFlags                   update_flags);
 
  template <bool level_dof_access>
  void
  reinit(
    const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell,
    const unsigned int face_no);
 
  template <bool level_dof_access>
  void
  reinit(
    const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell,
    const typename Triangulation<dim, spacedim>::face_iterator &          face);
 
  void
  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
         const unsigned int                                          face_no);
 
  /*
   * Reinitialize the gradients, Jacobi determinants, etc for the given face
   * on a given cell of type "iterator into a Triangulation object", and the
   * given finite element. Since iterators into a triangulation alone only
   * convey information about the geometry of a cell, but not about degrees of
   * freedom possibly associated with this cell, you will not be able to call
   * some functions of this class if they need information about degrees of
   * freedom. These functions are, above all, the
   * <tt>get_function_value/gradients/hessians/third_derivatives</tt>
   * functions. If you want to call these functions, you have to call the @p
   * reinit variants that take iterators into DoFHandler or other DoF handler
   * type objects.
   *
   * @note @p face must be one of @p cell's face iterators.
   */
  void
  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
         const typename Triangulation<dim, spacedim>::face_iterator &face);
 
  const FEFaceValues<dim, spacedim> &
  get_present_fe_values() const;
 
private:
  void
  initialize(const UpdateFlags update_flags);
 
  void
  do_reinit(const unsigned int face_no);
};
 
 
template <int dim, int spacedim = dim>
class FESubfaceValues : public FEFaceValuesBase<dim, spacedim>
{
public:
  static constexpr unsigned int dimension = dim;
 
  static constexpr unsigned int space_dimension = spacedim;
 
  static constexpr unsigned int integral_dimension = dim - 1;
 
  FESubfaceValues(const Mapping<dim, spacedim> &      mapping,
                  const FiniteElement<dim, spacedim> &fe,
                  const Quadrature<dim - 1> &         face_quadrature,
                  const UpdateFlags                   update_flags);
 
  FESubfaceValues(const Mapping<dim, spacedim> &      mapping,
                  const FiniteElement<dim, spacedim> &fe,
                  const hp::QCollection<dim - 1> &    face_quadrature,
                  const UpdateFlags                   update_flags);
 
  FESubfaceValues(const FiniteElement<dim, spacedim> &fe,
                  const Quadrature<dim - 1> &         face_quadrature,
                  const UpdateFlags                   update_flags);
 
  FESubfaceValues(const FiniteElement<dim, spacedim> &fe,
                  const hp::QCollection<dim - 1> &    face_quadrature,
                  const UpdateFlags                   update_flags);
 
  template <bool level_dof_access>
  void
  reinit(
    const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell,
    const unsigned int face_no,
    const unsigned int subface_no);
 
  template <bool level_dof_access>
  void
  reinit(
    const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell,
    const typename Triangulation<dim, spacedim>::face_iterator &          face,
    const typename Triangulation<dim, spacedim>::face_iterator &subface);
 
  void
  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
         const unsigned int                                          face_no,
         const unsigned int subface_no);
 
  void
  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
         const typename Triangulation<dim, spacedim>::face_iterator &face,
         const typename Triangulation<dim, spacedim>::face_iterator &subface);
 
  const FESubfaceValues<dim, spacedim> &
  get_present_fe_values() const;
 
  DeclException0(ExcReinitCalledWithBoundaryFace);
 
  DeclException0(ExcFaceHasNoSubfaces);
 
private:
  void
  initialize(const UpdateFlags update_flags);
 
  void
  do_reinit(const unsigned int face_no, const unsigned int subface_no);
};
 
 
#ifndef DOXYGEN
 
 
/*------------------------ Inline functions: namespace FEValuesViews --------*/
 
namespace FEValuesViews
{
  template <int dim, int spacedim>
  inline typename Scalar<dim, spacedim>::value_type
  Scalar<dim, spacedim>::value(const unsigned int shape_function,
                               const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(
      fe_values->update_flags & update_values,
      ((typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
        "update_values"))));
 
    // an adaptation of the FEValuesBase::shape_value_component function
    // except that here we know the component as fixed and we have
    // pre-computed and cached a bunch of information. See the comments there.
    if (shape_function_data[shape_function].is_nonzero_shape_function_component)
      return fe_values->finite_element_output.shape_values(
        shape_function_data[shape_function].row_index, q_point);
    else
      return 0;
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Scalar<dim, spacedim>::gradient_type
  Scalar<dim, spacedim>::gradient(const unsigned int shape_function,
                                  const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_gradients,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_gradients")));
 
    // an adaptation of the FEValuesBase::shape_grad_component
    // function except that here we know the component as fixed and we have
    // pre-computed and cached a bunch of information. See the comments there.
    if (shape_function_data[shape_function].is_nonzero_shape_function_component)
      return fe_values->finite_element_output
        .shape_gradients[shape_function_data[shape_function].row_index]
                        [q_point];
    else
      return gradient_type();
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Scalar<dim, spacedim>::hessian_type
  Scalar<dim, spacedim>::hessian(const unsigned int shape_function,
                                 const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_hessians,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_hessians")));
 
    // an adaptation of the FEValuesBase::shape_hessian_component
    // function except that here we know the component as fixed and we have
    // pre-computed and cached a bunch of information. See the comments there.
    if (shape_function_data[shape_function].is_nonzero_shape_function_component)
      return fe_values->finite_element_output
        .shape_hessians[shape_function_data[shape_function].row_index][q_point];
    else
      return hessian_type();
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Scalar<dim, spacedim>::third_derivative_type
  Scalar<dim, spacedim>::third_derivative(const unsigned int shape_function,
                                          const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_3rd_derivatives,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_3rd_derivatives")));
 
    // an adaptation of the FEValuesBase::shape_3rdderivative_component
    // function except that here we know the component as fixed and we have
    // pre-computed and cached a bunch of information. See the comments there.
    if (shape_function_data[shape_function].is_nonzero_shape_function_component)
      return fe_values->finite_element_output
        .shape_3rd_derivatives[shape_function_data[shape_function].row_index]
                              [q_point];
    else
      return third_derivative_type();
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Vector<dim, spacedim>::value_type
  Vector<dim, spacedim>::value(const unsigned int shape_function,
                               const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_values,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_values")));
 
    // same as for the scalar case except that we have one more index
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return value_type();
    else if (snc != -1)
      {
        value_type return_value;
        return_value[shape_function_data[shape_function]
                       .single_nonzero_component_index] =
          fe_values->finite_element_output.shape_values(snc, q_point);
        return return_value;
      }
    else
      {
        value_type return_value;
        for (unsigned int d = 0; d < dim; ++d)
          if (shape_function_data[shape_function]
                .is_nonzero_shape_function_component[d])
            return_value[d] = fe_values->finite_element_output.shape_values(
              shape_function_data[shape_function].row_index[d], q_point);
 
        return return_value;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Vector<dim, spacedim>::gradient_type
  Vector<dim, spacedim>::gradient(const unsigned int shape_function,
                                  const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_gradients,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_gradients")));
 
    // same as for the scalar case except that we have one more index
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return gradient_type();
    else if (snc != -1)
      {
        gradient_type return_value;
        return_value[shape_function_data[shape_function]
                       .single_nonzero_component_index] =
          fe_values->finite_element_output.shape_gradients[snc][q_point];
        return return_value;
      }
    else
      {
        gradient_type return_value;
        for (unsigned int d = 0; d < dim; ++d)
          if (shape_function_data[shape_function]
                .is_nonzero_shape_function_component[d])
            return_value[d] =
              fe_values->finite_element_output.shape_gradients
                [shape_function_data[shape_function].row_index[d]][q_point];
 
        return return_value;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Vector<dim, spacedim>::divergence_type
  Vector<dim, spacedim>::divergence(const unsigned int shape_function,
                                    const unsigned int q_point) const
  {
    // this function works like in the case above
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_gradients,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_gradients")));
 
    // same as for the scalar case except that we have one more index
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return divergence_type();
    else if (snc != -1)
      return fe_values->finite_element_output
        .shape_gradients[snc][q_point][shape_function_data[shape_function]
                                         .single_nonzero_component_index];
    else
      {
        divergence_type return_value = 0;
        for (unsigned int d = 0; d < dim; ++d)
          if (shape_function_data[shape_function]
                .is_nonzero_shape_function_component[d])
            return_value +=
              fe_values->finite_element_output.shape_gradients
                [shape_function_data[shape_function].row_index[d]][q_point][d];
 
        return return_value;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Vector<dim, spacedim>::curl_type
  Vector<dim, spacedim>::curl(const unsigned int shape_function,
                              const unsigned int q_point) const
  {
    // this function works like in the case above
 
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_gradients,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_gradients")));
    // same as for the scalar case except that we have one more index
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
 
    if (snc == -2)
      return curl_type();
 
    else
      switch (dim)
        {
          case 1:
            {
              Assert(false,
                     ExcMessage(
                       "Computing the curl in 1d is not a useful operation"));
              return curl_type();
            }
 
          case 2:
            {
              if (snc != -1)
                {
                  curl_type return_value;
 
                  // the single nonzero component can only be zero or one in 2d
                  if (shape_function_data[shape_function]
                        .single_nonzero_component_index == 0)
                    return_value[0] =
                      -1.0 * fe_values->finite_element_output
                               .shape_gradients[snc][q_point][1];
                  else
                    return_value[0] = fe_values->finite_element_output
                                        .shape_gradients[snc][q_point][0];
 
                  return return_value;
                }
 
              else
                {
                  curl_type return_value;
 
                  return_value[0] = 0.0;
 
                  if (shape_function_data[shape_function]
                        .is_nonzero_shape_function_component[0])
                    return_value[0] -=
                      fe_values->finite_element_output
                        .shape_gradients[shape_function_data[shape_function]
                                           .row_index[0]][q_point][1];
 
                  if (shape_function_data[shape_function]
                        .is_nonzero_shape_function_component[1])
                    return_value[0] +=
                      fe_values->finite_element_output
                        .shape_gradients[shape_function_data[shape_function]
                                           .row_index[1]][q_point][0];
 
                  return return_value;
                }
            }
 
          case 3:
            {
              if (snc != -1)
                {
                  curl_type return_value;
 
                  switch (shape_function_data[shape_function]
                            .single_nonzero_component_index)
                    {
                      case 0:
                        {
                          return_value[0] = 0;
                          return_value[1] = fe_values->finite_element_output
                                              .shape_gradients[snc][q_point][2];
                          return_value[2] =
                            -1.0 * fe_values->finite_element_output
                                     .shape_gradients[snc][q_point][1];
                          return return_value;
                        }
 
                      case 1:
                        {
                          return_value[0] =
                            -1.0 * fe_values->finite_element_output
                                     .shape_gradients[snc][q_point][2];
                          return_value[1] = 0;
                          return_value[2] = fe_values->finite_element_output
                                              .shape_gradients[snc][q_point][0];
                          return return_value;
                        }
 
                      default:
                        {
                          return_value[0] = fe_values->finite_element_output
                                              .shape_gradients[snc][q_point][1];
                          return_value[1] =
                            -1.0 * fe_values->finite_element_output
                                     .shape_gradients[snc][q_point][0];
                          return_value[2] = 0;
                          return return_value;
                        }
                    }
                }
 
              else
                {
                  curl_type return_value;
 
                  for (unsigned int i = 0; i < dim; ++i)
                    return_value[i] = 0.0;
 
                  if (shape_function_data[shape_function]
                        .is_nonzero_shape_function_component[0])
                    {
                      return_value[1] +=
                        fe_values->finite_element_output
                          .shape_gradients[shape_function_data[shape_function]
                                             .row_index[0]][q_point][2];
                      return_value[2] -=
                        fe_values->finite_element_output
                          .shape_gradients[shape_function_data[shape_function]
                                             .row_index[0]][q_point][1];
                    }
 
                  if (shape_function_data[shape_function]
                        .is_nonzero_shape_function_component[1])
                    {
                      return_value[0] -=
                        fe_values->finite_element_output
                          .shape_gradients[shape_function_data[shape_function]
                                             .row_index[1]][q_point][2];
                      return_value[2] +=
                        fe_values->finite_element_output
                          .shape_gradients[shape_function_data[shape_function]
                                             .row_index[1]][q_point][0];
                    }
 
                  if (shape_function_data[shape_function]
                        .is_nonzero_shape_function_component[2])
                    {
                      return_value[0] +=
                        fe_values->finite_element_output
                          .shape_gradients[shape_function_data[shape_function]
                                             .row_index[2]][q_point][1];
                      return_value[1] -=
                        fe_values->finite_element_output
                          .shape_gradients[shape_function_data[shape_function]
                                             .row_index[2]][q_point][0];
                    }
 
                  return return_value;
                }
            }
        }
    // should not end up here
    Assert(false, ExcInternalError());
    return curl_type();
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Vector<dim, spacedim>::hessian_type
  Vector<dim, spacedim>::hessian(const unsigned int shape_function,
                                 const unsigned int q_point) const
  {
    // this function works like in the case above
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_hessians,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_hessians")));
 
    // same as for the scalar case except that we have one more index
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return hessian_type();
    else if (snc != -1)
      {
        hessian_type return_value;
        return_value[shape_function_data[shape_function]
                       .single_nonzero_component_index] =
          fe_values->finite_element_output.shape_hessians[snc][q_point];
        return return_value;
      }
    else
      {
        hessian_type return_value;
        for (unsigned int d = 0; d < dim; ++d)
          if (shape_function_data[shape_function]
                .is_nonzero_shape_function_component[d])
            return_value[d] =
              fe_values->finite_element_output.shape_hessians
                [shape_function_data[shape_function].row_index[d]][q_point];
 
        return return_value;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Vector<dim, spacedim>::third_derivative_type
  Vector<dim, spacedim>::third_derivative(const unsigned int shape_function,
                                          const unsigned int q_point) const
  {
    // this function works like in the case above
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_3rd_derivatives,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_3rd_derivatives")));
 
    // same as for the scalar case except that we have one more index
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return third_derivative_type();
    else if (snc != -1)
      {
        third_derivative_type return_value;
        return_value[shape_function_data[shape_function]
                       .single_nonzero_component_index] =
          fe_values->finite_element_output.shape_3rd_derivatives[snc][q_point];
        return return_value;
      }
    else
      {
        third_derivative_type return_value;
        for (unsigned int d = 0; d < dim; ++d)
          if (shape_function_data[shape_function]
                .is_nonzero_shape_function_component[d])
            return_value[d] =
              fe_values->finite_element_output.shape_3rd_derivatives
                [shape_function_data[shape_function].row_index[d]][q_point];
 
        return return_value;
      }
  }
 
 
 
  namespace internal
  {
    inline ::SymmetricTensor<2, 1>
    symmetrize_single_row(const unsigned int n, const ::Tensor<1, 1> &t)
    {
      AssertIndexRange(n, 1);
      (void)n;
 
      return {{t[0]}};
    }
 
 
 
    inline ::SymmetricTensor<2, 2>
    symmetrize_single_row(const unsigned int n, const ::Tensor<1, 2> &t)
    {
      switch (n)
        {
          case 0:
            {
              return {{t[0], 0, t[1] / 2}};
            }
          case 1:
            {
              return {{0, t[1], t[0] / 2}};
            }
          default:
            {
              AssertIndexRange(n, 2);
              return {};
            }
        }
    }
 
 
 
    inline ::SymmetricTensor<2, 3>
    symmetrize_single_row(const unsigned int n, const ::Tensor<1, 3> &t)
    {
      switch (n)
        {
          case 0:
            {
              return {{t[0], 0, 0, t[1] / 2, t[2] / 2, 0}};
            }
          case 1:
            {
              return {{0, t[1], 0, t[0] / 2, 0, t[2] / 2}};
            }
          case 2:
            {
              return {{0, 0, t[2], 0, t[0] / 2, t[1] / 2}};
            }
          default:
            {
              AssertIndexRange(n, 3);
              return {};
            }
        }
    }
  } // namespace internal
 
 
 
  template <int dim, int spacedim>
  inline typename Vector<dim, spacedim>::symmetric_gradient_type
  Vector<dim, spacedim>::symmetric_gradient(const unsigned int shape_function,
                                            const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_gradients,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_gradients")));
 
    // same as for the scalar case except that we have one more index
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
    if (snc == -2)
      return symmetric_gradient_type();
    else if (snc != -1)
      return internal::symmetrize_single_row(
        shape_function_data[shape_function].single_nonzero_component_index,
        fe_values->finite_element_output.shape_gradients[snc][q_point]);
    else
      {
        gradient_type return_value;
        for (unsigned int d = 0; d < dim; ++d)
          if (shape_function_data[shape_function]
                .is_nonzero_shape_function_component[d])
            return_value[d] =
              fe_values->finite_element_output.shape_gradients
                [shape_function_data[shape_function].row_index[d]][q_point];
 
        return symmetrize(return_value);
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline typename SymmetricTensor<2, dim, spacedim>::value_type
  SymmetricTensor<2, dim, spacedim>::value(const unsigned int shape_function,
                                           const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_values,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_values")));
 
    // similar to the vector case where we have more then one index and we need
    // to convert between unrolled and component indexing for tensors
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
 
    if (snc == -2)
      {
        // shape function is zero for the selected components
        return value_type();
      }
    else if (snc != -1)
      {
        value_type         return_value;
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        return_value[value_type::unrolled_to_component_indices(comp)] =
          fe_values->finite_element_output.shape_values(snc, q_point);
        return return_value;
      }
    else
      {
        value_type return_value;
        for (unsigned int d = 0; d < value_type::n_independent_components; ++d)
          if (shape_function_data[shape_function]
                .is_nonzero_shape_function_component[d])
            return_value[value_type::unrolled_to_component_indices(d)] =
              fe_values->finite_element_output.shape_values(
                shape_function_data[shape_function].row_index[d], q_point);
        return return_value;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline typename SymmetricTensor<2, dim, spacedim>::divergence_type
  SymmetricTensor<2, dim, spacedim>::divergence(
    const unsigned int shape_function,
    const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_gradients,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_gradients")));
 
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
 
    if (snc == -2)
      {
        // shape function is zero for the selected components
        return divergence_type();
      }
    else if (snc != -1)
      {
        // we have a single non-zero component when the symmetric tensor is
        // represented in unrolled form. this implies we potentially have
        // two non-zero components when represented in component form!  we
        // will only have one non-zero entry if the non-zero component lies on
        // the diagonal of the tensor.
        //
        // the divergence of a second-order tensor is a first order tensor.
        //
        // assume the second-order tensor is A with components A_{ij}.  then
        // A_{ij} = A_{ji} and there is only one (if diagonal) or two non-zero
        // entries in the tensorial representation.  define the
        // divergence as:
        // b_i \dealcoloneq \dfrac{\partial phi_{ij}}{\partial x_j}.
        // (which is incidentally also
        // b_j \dealcoloneq \dfrac{\partial phi_{ij}}{\partial x_i}).
        // In both cases, a sum is implied.
        //
        // Now, we know the nonzero component in unrolled form: it is indicated
        // by 'snc'. we can figure out which tensor components belong to this:
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        const unsigned int ii =
          value_type::unrolled_to_component_indices(comp)[0];
        const unsigned int jj =
          value_type::unrolled_to_component_indices(comp)[1];
 
        // given the form of the divergence above, if ii=jj there is only a
        // single nonzero component of the full tensor and the gradient
        // equals
        // b_ii \dealcoloneq \dfrac{\partial phi_{ii,ii}}{\partial x_ii}.
        // all other entries of 'b' are zero
        //
        // on the other hand, if ii!=jj, then there are two nonzero entries in
        // the full tensor and
        // b_ii \dealcoloneq \dfrac{\partial phi_{ii,jj}}{\partial x_ii}.
        // b_jj \dealcoloneq \dfrac{\partial phi_{ii,jj}}{\partial x_jj}.
        // again, all other entries of 'b' are zero
        const ::Tensor<1, spacedim> &phi_grad =
          fe_values->finite_element_output.shape_gradients[snc][q_point];
 
        divergence_type return_value;
        return_value[ii] = phi_grad[jj];
 
        if (ii != jj)
          return_value[jj] = phi_grad[ii];
 
        return return_value;
      }
    else
      {
        Assert(false, ExcNotImplemented());
        divergence_type return_value;
        return return_value;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Tensor<2, dim, spacedim>::value_type
  Tensor<2, dim, spacedim>::value(const unsigned int shape_function,
                                  const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_values,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_values")));
 
    // similar to the vector case where we have more then one index and we need
    // to convert between unrolled and component indexing for tensors
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
 
    if (snc == -2)
      {
        // shape function is zero for the selected components
        return value_type();
      }
    else if (snc != -1)
      {
        value_type         return_value;
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        const TableIndices<2> indices =
          ::Tensor<2, spacedim>::unrolled_to_component_indices(comp);
        return_value[indices] =
          fe_values->finite_element_output.shape_values(snc, q_point);
        return return_value;
      }
    else
      {
        value_type return_value;
        for (unsigned int d = 0; d < dim * dim; ++d)
          if (shape_function_data[shape_function]
                .is_nonzero_shape_function_component[d])
            {
              const TableIndices<2> indices =
                ::Tensor<2, spacedim>::unrolled_to_component_indices(d);
              return_value[indices] =
                fe_values->finite_element_output.shape_values(
                  shape_function_data[shape_function].row_index[d], q_point);
            }
        return return_value;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Tensor<2, dim, spacedim>::divergence_type
  Tensor<2, dim, spacedim>::divergence(const unsigned int shape_function,
                                       const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_gradients,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_gradients")));
 
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
 
    if (snc == -2)
      {
        // shape function is zero for the selected components
        return divergence_type();
      }
    else if (snc != -1)
      {
        // we have a single non-zero component when the tensor is
        // represented in unrolled form.
        //
        // the divergence of a second-order tensor is a first order tensor.
        //
        // assume the second-order tensor is A with components A_{ij},
        // then divergence is d_i := \frac{\partial A_{ij}}{\partial x_j}
        //
        // Now, we know the nonzero component in unrolled form: it is indicated
        // by 'snc'. we can figure out which tensor components belong to this:
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        const TableIndices<2> indices =
          ::Tensor<2, spacedim>::unrolled_to_component_indices(comp);
        const unsigned int ii = indices[0];
        const unsigned int jj = indices[1];
 
        const ::Tensor<1, spacedim> &phi_grad =
          fe_values->finite_element_output.shape_gradients[snc][q_point];
 
        divergence_type return_value;
        // note that we contract \nabla from the right
        return_value[ii] = phi_grad[jj];
 
        return return_value;
      }
    else
      {
        Assert(false, ExcNotImplemented());
        divergence_type return_value;
        return return_value;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline typename Tensor<2, dim, spacedim>::gradient_type
  Tensor<2, dim, spacedim>::gradient(const unsigned int shape_function,
                                     const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
    Assert(fe_values->update_flags & update_gradients,
           (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
             "update_gradients")));
 
    const int snc =
      shape_function_data[shape_function].single_nonzero_component;
 
    if (snc == -2)
      {
        // shape function is zero for the selected components
        return gradient_type();
      }
    else if (snc != -1)
      {
        // we have a single non-zero component when the tensor is
        // represented in unrolled form.
        //
        // the gradient of a second-order tensor is a third order tensor.
        //
        // assume the second-order tensor is A with components A_{ij},
        // then gradient is B_{ijk} := \frac{\partial A_{ij}}{\partial x_k}
        //
        // Now, we know the nonzero component in unrolled form: it is indicated
        // by 'snc'. we can figure out which tensor components belong to this:
        const unsigned int comp =
          shape_function_data[shape_function].single_nonzero_component_index;
        const TableIndices<2> indices =
          ::Tensor<2, spacedim>::unrolled_to_component_indices(comp);
        const unsigned int ii = indices[0];
        const unsigned int jj = indices[1];
 
        const ::Tensor<1, spacedim> &phi_grad =
          fe_values->finite_element_output.shape_gradients[snc][q_point];
 
        gradient_type return_value;
        return_value[ii][jj] = phi_grad;
 
        return return_value;
      }
    else
      {
        Assert(false, ExcNotImplemented());
        gradient_type return_value;
        return return_value;
      }
  }
 
} // namespace FEValuesViews
 
 
 
/*---------------------- Inline functions: FEValuesBase ---------------------*/
 
 
 
template <int dim, int spacedim>
template <bool lda>
inline FEValuesBase<dim, spacedim>::CellIteratorContainer::
  CellIteratorContainer(
    const TriaIterator<DoFCellAccessor<dim, spacedim, lda>> &cell)
  : initialized(true)
  , cell(cell)
  , dof_handler(&cell->get_dof_handler())
  , level_dof_access(lda)
{}
 
 
 
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];
}
 
 
 
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];
}
 
 
 
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];
}
 
 
 
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];
}
 
 
 
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];
}
 
 
 
/*--------------------- Inline functions: FEValues --------------------------*/
 
 
template <int dim, int spacedim>
inline const Quadrature<dim> &
FEValues<dim, spacedim>::get_quadrature() const
{
  return quadrature;
}
 
 
 
template <int dim, int spacedim>
inline const FEValues<dim, spacedim> &
FEValues<dim, spacedim>::get_present_fe_values() const
{
  return *this;
}
 
 
/*---------------------- Inline functions: FEFaceValuesBase -----------------*/
 
 
template <int dim, int spacedim>
inline unsigned int
FEFaceValuesBase<dim, spacedim>::get_face_number() const
{
  return present_face_no;
}
 
 
template <int dim, int spacedim>
inline unsigned int
FEFaceValuesBase<dim, spacedim>::get_face_index() const
{
  return present_face_index;
}
 
 
/*----------------------- Inline functions: FE*FaceValues -------------------*/
 
template <int dim, int spacedim>
inline const Quadrature<dim - 1> &
FEFaceValuesBase<dim, spacedim>::get_quadrature() const
{
  return quadrature[quadrature.size() == 1 ? 0 : present_face_no];
}
 
 
 
template <int dim, int spacedim>
inline const FEFaceValues<dim, spacedim> &
FEFaceValues<dim, spacedim>::get_present_fe_values() const
{
  return *this;
}
 
 
 
template <int dim, int spacedim>
inline const FESubfaceValues<dim, spacedim> &
FESubfaceValues<dim, spacedim>::get_present_fe_values() const
{
  return *this;
}
 
 
 
template <int dim, int spacedim>
inline const Tensor<1, spacedim> &
FEFaceValuesBase<dim, spacedim>::boundary_form(const unsigned int q_point) const
{
  AssertIndexRange(q_point, this->mapping_output.boundary_forms.size());
  Assert(this->update_flags & update_boundary_forms,
         (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
           "update_boundary_forms")));
 
  return this->mapping_output.boundary_forms[q_point];
}
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif