fe_evaluation.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2012 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
 
#ifndef dealii_matrix_free_fe_evaluation_h
#define dealii_matrix_free_fe_evaluation_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/array_view.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/observer_pointer.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/vectorization.h>
 
#include <deal.II/fe/fe_values_views.h>
 
#include <deal.II/lac/vector_operation.h>
 
#include <deal.II/matrix_free/evaluation_flags.h>
#include <deal.II/matrix_free/evaluation_kernels.h>
#include <deal.II/matrix_free/evaluation_kernels_face.h>
#include <deal.II/matrix_free/evaluation_template_factory.h>
#include <deal.II/matrix_free/fe_evaluation_data.h>
#include <deal.II/matrix_free/hanging_nodes_internal.h>
#include <deal.II/matrix_free/mapping_data_on_the_fly.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/tensor_product_kernels.h>
#include <deal.II/matrix_free/type_traits.h>
#include <deal.II/matrix_free/vector_access_internal.h>
 
#include <type_traits>
 
 
DEAL_II_NAMESPACE_OPEN
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
class FEEvaluationBase
  : public FEEvaluationData<dim, VectorizedArrayType, is_face>
{
public:
  using number_type = Number;
  using value_type =
    std::conditional_t<n_components_ == 1,
                       VectorizedArrayType,
                       Tensor<1, n_components_, VectorizedArrayType>>;
  using gradient_type = std::conditional_t<
    n_components_ == 1,
    Tensor<1, dim, VectorizedArrayType>,
    std::conditional_t<
      n_components_ == dim,
      Tensor<2, dim, VectorizedArrayType>,
      Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>>>;
  using curl_type    = internal::CurlType<dim, VectorizedArrayType>;
  using hessian_type = std::conditional_t<
    n_components_ == 1,
    Tensor<2, dim, VectorizedArrayType>,
    std::conditional_t<
      n_components_ == dim,
      Tensor<3, dim, VectorizedArrayType>,
      Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>>>;
  static constexpr unsigned int dimension    = dim;
  static constexpr unsigned int n_components = n_components_;
  static constexpr unsigned int n_lanes      = VectorizedArrayType::size();
 
  template <typename VectorType>
  void
  read_dof_values(
    const VectorType           &src,
    const unsigned int          first_index = 0,
    const std::bitset<n_lanes> &mask        = std::bitset<n_lanes>().flip());
 
  template <typename VectorType>
  void
  read_dof_values_plain(
    const VectorType           &src,
    const unsigned int          first_index = 0,
    const std::bitset<n_lanes> &mask        = std::bitset<n_lanes>().flip());
 
  template <typename VectorType>
  void
  distribute_local_to_global(
    VectorType                 &dst,
    const unsigned int          first_index = 0,
    const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
 
  template <typename VectorType>
  void
  set_dof_values(
    VectorType                 &dst,
    const unsigned int          first_index = 0,
    const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
 
  template <typename VectorType>
  void
  set_dof_values_plain(
    VectorType                 &dst,
    const unsigned int          first_index = 0,
    const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
 
  value_type
  get_dof_value(const unsigned int dof) const;
 
  void
  submit_dof_value(const value_type val_in, const unsigned int dof);
 
  value_type
  get_value(const unsigned int q_point) const;
 
  void
  submit_value(const value_type val_in, const unsigned int q_point);
 
  template <int n_components_local = n_components,
            typename = std::enable_if_t<n_components == n_components_local>>
  void
  submit_value(const Tensor<1, 1, VectorizedArrayType> val_in,
               const unsigned int                      q_point);
 
  gradient_type
  get_gradient(const unsigned int q_point) const;
 
  value_type
  get_normal_derivative(const unsigned int q_point) const;
 
  void
  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
 
  template <int dim_ = dim,
            typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
  void
  submit_gradient(const Tensor<2, 1, VectorizedArrayType> val_in,
                  const unsigned int                      q_point);
 
  void
  submit_normal_derivative(const value_type   grad_in,
                           const unsigned int q_point);
 
  hessian_type
  get_hessian(const unsigned int q_point) const;
 
  gradient_type
  get_hessian_diagonal(const unsigned int q_point) const;
 
  value_type
  get_laplacian(const unsigned int q_point) const;
 
  value_type
  get_normal_hessian(const unsigned int q_point) const;
 
  void
  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
 
  void
  submit_normal_hessian(const value_type   normal_hessian_in,
                        const unsigned int q_point);
 
  template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
  VectorizedArrayType
  get_divergence(const unsigned int q_point) const;
 
  template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
  void
  submit_divergence(const VectorizedArrayType div_in,
                    const unsigned int        q_point);
 
  template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
  SymmetricTensor<2, dim, VectorizedArrayType>
  get_symmetric_gradient(const unsigned int q_point) const;
 
  template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
  void
  submit_symmetric_gradient(
    const SymmetricTensor<2, dim, VectorizedArrayType> grad_in,
    const unsigned int                                 q_point);
 
  template <int dim_ = dim,
            typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
  curl_type
  get_curl(const unsigned int q_point) const;
 
  template <int dim_ = dim,
            typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
  void
  submit_curl(const curl_type curl_in, const unsigned int q_point);
 
  value_type
  integrate_value() const;
 
  const MatrixFree<dim, Number, VectorizedArrayType> &
  get_matrix_free() const;
 
protected:
  FEEvaluationBase(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                  dof_handler_index,
    const unsigned int first_selected_component,
    const unsigned int quadrature_index,
    const unsigned int fe_degree,
    const unsigned int n_q_points,
    const bool         is_interior_face,
    const unsigned int active_fe_index,
    const unsigned int active_quad_index,
    const unsigned int face_type);
 
  FEEvaluationBase(
    const Mapping<dim>       &mapping,
    const FiniteElement<dim> &fe,
    const Quadrature<1>      &quadrature,
    const UpdateFlags         update_flags,
    const unsigned int        first_selected_component,
    const FEEvaluationData<dim, VectorizedArrayType, is_face> *other);
 
  FEEvaluationBase(const FEEvaluationBase &other);
 
  FEEvaluationBase &
  operator=(const FEEvaluationBase &other);
 
  ~FEEvaluationBase();
 
  template <typename VectorType, typename VectorOperation>
  void
  read_write_operation(
    const VectorOperation                         &operation,
    const std::array<VectorType *, n_components_> &vectors,
    const std::array<
      const std::vector<ArrayView<const typename VectorType::value_type>> *,
      n_components_>           &vectors_sm,
    const std::bitset<n_lanes> &mask,
    const bool                  apply_constraints = true) const;
 
  template <typename VectorType, typename VectorOperation>
  void
  read_write_operation_contiguous(
    const VectorOperation                         &operation,
    const std::array<VectorType *, n_components_> &vectors,
    const std::array<
      const std::vector<ArrayView<const typename VectorType::value_type>> *,
      n_components_>           &vectors_sm,
    const std::bitset<n_lanes> &mask) const;
 
  template <typename VectorType, typename VectorOperation>
  void
  read_write_operation_global(
    const VectorOperation                         &operation,
    const std::array<VectorType *, n_components_> &vectors) const;
 
  void
  apply_hanging_node_constraints(const bool transpose) const;
 
  AlignedVector<VectorizedArrayType> *scratch_data_array;
 
  const MatrixFree<dim, Number, VectorizedArrayType> *matrix_free;
 
  mutable std::vector<types::global_dof_index> local_dof_indices;
};
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
class FEEvaluation : public FEEvaluationBase<dim,
                                             n_components_,
                                             Number,
                                             false,
                                             VectorizedArrayType>
{
  static_assert(
    std::is_same_v<Number, typename VectorizedArrayType::value_type>,
    "Type of Number and of VectorizedArrayType do not match.");
 
public:
  using BaseClass =
    FEEvaluationBase<dim, n_components_, Number, false, VectorizedArrayType>;
 
  using number_type = Number;
 
  using value_type = typename BaseClass::value_type;
 
  using gradient_type = typename BaseClass::gradient_type;
 
  static constexpr unsigned int dimension = dim;
 
  static constexpr unsigned int n_components = n_components_;
 
  static constexpr unsigned int n_lanes = VectorizedArrayType::size();
 
  static constexpr unsigned int static_n_q_points =
    Utilities::pow(n_q_points_1d, dim);
 
  static constexpr unsigned int static_dofs_per_component =
    Utilities::pow(fe_degree + 1, dim);
 
  static constexpr unsigned int tensor_dofs_per_cell =
    static_dofs_per_component * n_components;
 
  static constexpr unsigned int static_dofs_per_cell =
    static_dofs_per_component * n_components;
 
  FEEvaluation(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                  dof_handler_index = 0,
    const unsigned int                                  quadrature_index  = 0,
    const unsigned int first_selected_component                           = 0,
    const unsigned int active_fe_index   = numbers::invalid_unsigned_int,
    const unsigned int active_quad_index = numbers::invalid_unsigned_int);
 
  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
               const std::pair<unsigned int, unsigned int>        &range,
               const unsigned int dof_handler_index        = 0,
               const unsigned int quadrature_index         = 0,
               const unsigned int first_selected_component = 0);
 
  FEEvaluation(const Mapping<dim>       &mapping,
               const FiniteElement<dim> &fe,
               const Quadrature<1>      &quadrature,
               const UpdateFlags         update_flags,
               const unsigned int        first_selected_component = 0);
 
  FEEvaluation(const FiniteElement<dim> &fe,
               const Quadrature<1>      &quadrature,
               const UpdateFlags         update_flags,
               const unsigned int        first_selected_component = 0);
 
  FEEvaluation(const FiniteElement<dim>                                &fe,
               const FEEvaluationData<dim, VectorizedArrayType, false> &other,
               const unsigned int first_selected_component = 0);
 
  FEEvaluation(const FEEvaluation &other);
 
  FEEvaluation &
  operator=(const FEEvaluation &other);
 
  void
  reinit(const unsigned int cell_batch_index);
 
  void
  reinit(const std::array<unsigned int, n_lanes> &cell_ids);
 
  template <bool level_dof_access>
  void
  reinit(const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell);
 
  void
  reinit(const typename Triangulation<dim>::cell_iterator &cell);
 
  static bool
  fast_evaluation_supported(const unsigned int given_degree,
                            const unsigned int given_n_q_points_1d);
 
  void
  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  void
  evaluate(const VectorizedArrayType             *values_array,
           const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  template <typename VectorType>
  void
  gather_evaluate(const VectorType                      &input_vector,
                  const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  void
  integrate(const EvaluationFlags::EvaluationFlags integration_flag);
 
  void
  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
            VectorizedArrayType                   *values_array,
            const bool                             sum_into_values = false);
 
  template <typename VectorType>
  void
  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
                    VectorType                            &output_vector);
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  dof_indices() const;
 
  const unsigned int dofs_per_component;
 
  const unsigned int dofs_per_cell;
 
  const unsigned int n_q_points;
 
private:
  void
  check_template_arguments(const unsigned int fe_no,
                           const unsigned int first_selected_component);
};
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d            = fe_degree + 1,
          int n_components_            = 1,
          typename Number              = double,
          typename VectorizedArrayType = VectorizedArray<Number>>
class FEFaceEvaluation : public FEEvaluationBase<dim,
                                                 n_components_,
                                                 Number,
                                                 true,
                                                 VectorizedArrayType>
{
  static_assert(
    std::is_same_v<Number, typename VectorizedArrayType::value_type>,
    "Type of Number and of VectorizedArrayType do not match.");
 
public:
  using BaseClass =
    FEEvaluationBase<dim, n_components_, Number, true, VectorizedArrayType>;
 
  using number_type = Number;
 
  using value_type = typename BaseClass::value_type;
 
  using gradient_type = typename BaseClass::gradient_type;
 
  static constexpr unsigned int dimension = dim;
 
  static constexpr unsigned int n_components = n_components_;
 
  static constexpr unsigned int n_lanes = VectorizedArrayType::size();
 
  static constexpr unsigned int static_n_q_points =
    Utilities::pow(n_q_points_1d, dim - 1);
 
  static constexpr unsigned int static_n_q_points_cell =
    Utilities::pow(n_q_points_1d, dim);
 
  static constexpr unsigned int static_dofs_per_component =
    Utilities::pow(fe_degree + 1, dim);
 
  static constexpr unsigned int tensor_dofs_per_cell =
    static_dofs_per_component * n_components;
 
  static constexpr unsigned int static_dofs_per_cell =
    static_dofs_per_component * n_components;
 
  FEFaceEvaluation(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const bool                                          is_interior_face = true,
    const unsigned int                                  dof_handler_index = 0,
    const unsigned int                                  quadrature_index  = 0,
    const unsigned int first_selected_component                           = 0,
    const unsigned int active_fe_index   = numbers::invalid_unsigned_int,
    const unsigned int active_quad_index = numbers::invalid_unsigned_int,
    const unsigned int face_type         = numbers::invalid_unsigned_int);
 
  FEFaceEvaluation(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const std::pair<unsigned int, unsigned int>        &range,
    const bool                                          is_interior_face = true,
    const unsigned int                                  dof_handler_index = 0,
    const unsigned int                                  quadrature_index  = 0,
    const unsigned int first_selected_component                           = 0);
 
  void
  reinit(const unsigned int face_batch_number);
 
  void
  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
 
  static bool
  fast_evaluation_supported(const unsigned int given_degree,
                            const unsigned int given_n_q_points_1d);
 
  void
  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  void
  evaluate(const VectorizedArrayType             *values_array,
           const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  void
  project_to_face(const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  void
  project_to_face(const VectorizedArrayType             *values_array,
                  const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  void
  evaluate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  template <typename VectorType>
  void
  gather_evaluate(const VectorType                      &input_vector,
                  const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  void
  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
            const bool                             sum_into_values = false);
 
  void
  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
            VectorizedArrayType                   *values_array,
            const bool                             sum_into_values = false);
 
  void
  integrate_in_face(const EvaluationFlags::EvaluationFlags integration_flag);
 
  void
  collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag,
                    const bool sum_into_values = false);
 
  void
  collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag,
                    VectorizedArrayType                   *values_array,
                    const bool sum_into_values = false);
 
  template <typename VectorType>
  void
  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
                    VectorType                            &output_vector);
 
  template <typename VectorType>
  void
  integrate_scatter(const bool  integrate_values,
                    const bool  integrate_gradients,
                    VectorType &output_vector);
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  dof_indices() const;
 
  bool
  at_boundary() const;
 
  types::boundary_id
  boundary_id() const;
 
  unsigned int
  get_dofs_per_component_projected_to_face();
 
  unsigned int
  get_dofs_projected_to_face();
 
  const unsigned int dofs_per_component;
 
  const unsigned int dofs_per_cell;
 
  const unsigned int n_q_points;
};
 
 
 
namespace internal
{
  namespace MatrixFreeFunctions
  {
    // a helper function to compute the number of DoFs of a DGP element at
    // compile time, depending on the degree
    template <int dim, int degree>
    struct DGP_dofs_per_component
    {
      // this division is always without remainder
      static constexpr unsigned int value =
        (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
    };
 
    // base specialization: 1d elements have 'degree+1' degrees of freedom
    template <int degree>
    struct DGP_dofs_per_component<1, degree>
    {
      static constexpr unsigned int value = degree + 1;
    };
  } // namespace MatrixFreeFunctions
} // namespace internal
 
 
/*----------------------- Inline functions ----------------------------------*/
 
#ifndef DOXYGEN
 
 
namespace internal
{
  // Extract all internal data pointers and indices in a single function that
  // get passed on to the constructor of FEEvaluationData, avoiding to look
  // things up multiple times
  template <bool is_face,
            int  dim,
            typename Number,
            typename VectorizedArrayType>
  inline typename FEEvaluationData<dim, VectorizedArrayType, is_face>::
    InitializationData
    extract_initialization_data(
      const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
      const unsigned int                                  dof_handler_index,
      const unsigned int first_selected_component,
      const unsigned int quadrature_index,
      const unsigned int fe_degree,
      const unsigned int n_q_points,
      const unsigned int active_fe_index_given,
      const unsigned int active_quad_index_given,
      const unsigned int face_type)
  {
    typename FEEvaluationData<dim, VectorizedArrayType, is_face>::
      InitializationData init_data;
 
    init_data.dof_info = &matrix_free.get_dof_info(dof_handler_index);
    init_data.mapping_data =
      &internal::MatrixFreeFunctions::
        MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
          matrix_free.get_mapping_info(), quadrature_index);
 
    init_data.active_fe_index =
      fe_degree != numbers::invalid_unsigned_int ?
        init_data.dof_info->fe_index_from_degree(first_selected_component,
                                                 fe_degree) :
        (active_fe_index_given != numbers::invalid_unsigned_int ?
           active_fe_index_given :
           0);
 
    init_data.active_quad_index =
      fe_degree == numbers::invalid_unsigned_int ?
        (active_quad_index_given != numbers::invalid_unsigned_int ?
           active_quad_index_given :
           std::min<unsigned int>(
             init_data.active_fe_index,
             init_data.mapping_data->descriptor.size() /
                 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
               1)) :
        init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
 
    init_data.shape_info = &matrix_free.get_shape_info(
      dof_handler_index,
      quadrature_index,
      init_data.dof_info->component_to_base_index[first_selected_component],
      init_data.active_fe_index,
      init_data.active_quad_index);
    init_data.descriptor =
      &init_data.mapping_data->descriptor
         [is_face ?
            (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
             (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
            init_data.active_quad_index];
 
    return init_data;
  }
} // namespace internal
 
 
 
/*----------------------- FEEvaluationBase ----------------------------------*/
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline FEEvaluationBase<dim,
                        n_components_,
                        Number,
                        is_face,
                        VectorizedArrayType>::
  FEEvaluationBase(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                  dof_handler_index,
    const unsigned int first_selected_component,
    const unsigned int quadrature_index,
    const unsigned int fe_degree,
    const unsigned int n_q_points,
    const bool         is_interior_face,
    const unsigned int active_fe_index,
    const unsigned int active_quad_index,
    const unsigned int face_type)
  : FEEvaluationData<dim, VectorizedArrayType, is_face>(
      internal::extract_initialization_data<is_face>(matrix_free,
                                                     dof_handler_index,
                                                     first_selected_component,
                                                     quadrature_index,
                                                     fe_degree,
                                                     n_q_points,
                                                     active_fe_index,
                                                     active_quad_index,
                                                     face_type),
      is_interior_face,
      quadrature_index,
      first_selected_component)
  , scratch_data_array(matrix_free.acquire_scratch_data())
  , matrix_free(&matrix_free)
{
  this->set_data_pointers(scratch_data_array, n_components_);
  Assert(
    this->dof_info->start_components.back() == 1 ||
      static_cast<int>(n_components_) <=
        static_cast<int>(
          this->dof_info->start_components
            [this->dof_info->component_to_base_index[first_selected_component] +
             1]) -
          first_selected_component,
    ExcMessage(
      "You tried to construct a vector-valued evaluator with " +
      std::to_string(n_components) +
      " components. However, "
      "the current base element has only " +
      std::to_string(
        this->dof_info->start_components
          [this->dof_info->component_to_base_index[first_selected_component] +
           1] -
        first_selected_component) +
      " components left when starting from local element index " +
      std::to_string(
        first_selected_component -
        this->dof_info->start_components
          [this->dof_info->component_to_base_index[first_selected_component]]) +
      " (global index " + std::to_string(first_selected_component) + ")"));
 
  // do not check for correct dimensions of data fields here, should be done
  // in derived classes
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline FEEvaluationBase<dim,
                        n_components_,
                        Number,
                        is_face,
                        VectorizedArrayType>::
  FEEvaluationBase(
    const Mapping<dim>       &mapping,
    const FiniteElement<dim> &fe,
    const Quadrature<1>      &quadrature,
    const UpdateFlags         update_flags,
    const unsigned int        first_selected_component,
    const FEEvaluationData<dim, VectorizedArrayType, is_face> *other)
  : FEEvaluationData<dim, VectorizedArrayType, is_face>(
      other != nullptr &&
          other->mapped_geometry->get_quadrature() == quadrature ?
        other->mapped_geometry :
        std::make_shared<internal::MatrixFreeFunctions::
                           MappingDataOnTheFly<dim, VectorizedArrayType>>(
          mapping,
          quadrature,
          update_flags),
      n_components_,
      first_selected_component)
  , scratch_data_array(new AlignedVector<VectorizedArrayType>())
  , matrix_free(nullptr)
{
  const unsigned int base_element_number =
    fe.component_to_base_index(first_selected_component).first;
  Assert(fe.element_multiplicity(base_element_number) == 1 ||
           fe.element_multiplicity(base_element_number) -
               first_selected_component >=
             n_components_,
         ExcMessage("The underlying element must at least contain as many "
                    "components as requested by this class"));
  (void)base_element_number;
 
  Assert(this->data == nullptr, ExcInternalError());
  this->data = new internal::MatrixFreeFunctions::ShapeInfo<Number>(
    Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
    fe,
    fe.component_to_base_index(first_selected_component).first);
 
  this->set_data_pointers(scratch_data_array, n_components_);
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline FEEvaluationBase<dim,
                        n_components_,
                        Number,
                        is_face,
                        VectorizedArrayType>::
  FEEvaluationBase(const FEEvaluationBase<dim,
                                          n_components_,
                                          Number,
                                          is_face,
                                          VectorizedArrayType> &other)
  : FEEvaluationData<dim, VectorizedArrayType, is_face>(other)
  , scratch_data_array(other.matrix_free == nullptr ?
                         new AlignedVector<VectorizedArrayType>() :
                         other.matrix_free->acquire_scratch_data())
  , matrix_free(other.matrix_free)
{
  if (other.matrix_free == nullptr)
    {
      Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
      this->data =
        new internal::MatrixFreeFunctions::ShapeInfo<Number>(*other.data);
 
      // Create deep copy of mapped geometry for use in parallel
      this->mapped_geometry =
        std::make_shared<internal::MatrixFreeFunctions::
                           MappingDataOnTheFly<dim, VectorizedArrayType>>(
          other.mapped_geometry->get_fe_values().get_mapping(),
          other.mapped_geometry->get_quadrature(),
          other.mapped_geometry->get_fe_values().get_update_flags());
 
      if constexpr (is_face == false)
        this->mapping_data = &this->mapped_geometry->get_data_storage();
      else
        Assert(false,
               ExcNotImplemented("On-the-fly geometry evaluation with "
                                 "face evaluators is not currently "
                                 "implemented!"));
 
      this->cell = 0;
 
      this->jacobian =
        this->mapped_geometry->get_data_storage().jacobians[0].begin();
      this->J_value =
        this->mapped_geometry->get_data_storage().JxW_values.begin();
      this->jacobian_gradients =
        this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
      this->jacobian_gradients_non_inverse =
        this->mapped_geometry->get_data_storage()
          .jacobian_gradients_non_inverse[0]
          .begin();
      this->quadrature_points =
        this->mapped_geometry->get_data_storage().quadrature_points.begin();
    }
 
  this->set_data_pointers(scratch_data_array, n_components_);
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline FEEvaluationBase<dim,
                        n_components_,
                        Number,
                        is_face,
                        VectorizedArrayType> &
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
operator=(const FEEvaluationBase<dim,
                                 n_components_,
                                 Number,
                                 is_face,
                                 VectorizedArrayType> &other)
{
  // release old memory
  if (matrix_free == nullptr)
    {
      delete this->data;
      delete scratch_data_array;
    }
  else
    {
      matrix_free->release_scratch_data(scratch_data_array);
    }
 
  this->FEEvaluationData<dim, VectorizedArrayType, is_face>::operator=(other);
 
  matrix_free = other.matrix_free;
 
  if (other.matrix_free == nullptr)
    {
      Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
      this->data =
        new internal::MatrixFreeFunctions::ShapeInfo<Number>(*other.data);
      scratch_data_array = new AlignedVector<VectorizedArrayType>();
 
      // Create deep copy of mapped geometry for use in parallel
      this->mapped_geometry =
        std::make_shared<internal::MatrixFreeFunctions::
                           MappingDataOnTheFly<dim, VectorizedArrayType>>(
          other.mapped_geometry->get_fe_values().get_mapping(),
          other.mapped_geometry->get_quadrature(),
          other.mapped_geometry->get_fe_values().get_update_flags());
 
      if constexpr (is_face == false)
        this->mapping_data = &this->mapped_geometry->get_data_storage();
      else
        Assert(false,
               ExcNotImplemented("On-the-fly geometry evaluation with "
                                 "face evaluators is not currently "
                                 "implemented!"));
      this->cell = 0;
 
      this->jacobian =
        this->mapped_geometry->get_data_storage().jacobians[0].begin();
      this->J_value =
        this->mapped_geometry->get_data_storage().JxW_values.begin();
      this->jacobian_gradients =
        this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
      this->jacobian_gradients_non_inverse =
        this->mapped_geometry->get_data_storage()
          .jacobian_gradients_non_inverse[0]
          .begin();
      this->quadrature_points =
        this->mapped_geometry->get_data_storage().quadrature_points.begin();
    }
  else
    {
      scratch_data_array = matrix_free->acquire_scratch_data();
    }
 
  this->set_data_pointers(scratch_data_array, n_components_);
 
  return *this;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline FEEvaluationBase<dim,
                        n_components_,
                        Number,
                        is_face,
                        VectorizedArrayType>::~FEEvaluationBase()
{
  if (matrix_free != nullptr)
    {
      try
        {
          matrix_free->release_scratch_data(scratch_data_array);
        }
      catch (...)
        {}
    }
  else
    {
      delete scratch_data_array;
      delete this->data;
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline const MatrixFree<dim, Number, VectorizedArrayType> &
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_matrix_free() const
{
  Assert(matrix_free != nullptr,
         ExcMessage(
           "FEEvaluation was not initialized with a MatrixFree object!"));
  return *matrix_free;
}
 
 
 
namespace internal
{
  // given a block vector return the underlying vector type
  // including constness (specified by bool)
  template <typename VectorType, bool>
  struct ConstBlockVectorSelector;
 
  template <typename VectorType>
  struct ConstBlockVectorSelector<VectorType, true>
  {
    using BaseVectorType = const typename VectorType::BlockType;
  };
 
  template <typename VectorType>
  struct ConstBlockVectorSelector<VectorType, false>
  {
    using BaseVectorType = typename VectorType::BlockType;
  };
 
  // allows to select between block vectors and non-block vectors, which
  // allows to use a unified interface for extracting blocks on block vectors
  // and doing nothing on usual vectors
  template <typename VectorType, bool>
  struct BlockVectorSelector;
 
  template <typename VectorType>
  struct BlockVectorSelector<VectorType, true>
  {
    using BaseVectorType = typename ConstBlockVectorSelector<
      VectorType,
      std::is_const_v<VectorType>>::BaseVectorType;
 
    static BaseVectorType *
    get_vector_component(VectorType &vec, const unsigned int component)
    {
      AssertIndexRange(component, vec.n_blocks());
      return &vec.block(component);
    }
  };
 
  template <typename VectorType>
  struct BlockVectorSelector<VectorType, false>
  {
    using BaseVectorType = VectorType;
 
    static BaseVectorType *
    get_vector_component(VectorType &vec, const unsigned int component)
    {
      // FEEvaluation allows to combine several vectors from a scalar
      // FiniteElement into a "vector-valued" FEEvaluation object with
      // multiple components. These components can be extracted with the other
      // get_vector_component functions. If we do not get a vector of vectors
      // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
      // must make sure that we do not duplicate the components in input
      // and/or duplicate the resulting integrals. In such a case, we should
      // only get the zeroth component in the vector contained set nullptr for
      // the others which allows us to catch unintended use in
      // read_write_operation.
      if (component == 0)
        return &vec;
      else
        return nullptr;
    }
  };
 
  template <typename VectorType>
  struct BlockVectorSelector<std::vector<VectorType>, false>
  {
    using BaseVectorType = VectorType;
 
    static BaseVectorType *
    get_vector_component(std::vector<VectorType> &vec,
                         const unsigned int       component)
    {
      AssertIndexRange(component, vec.size());
      return &vec[component];
    }
  };
 
  template <typename VectorType>
  struct BlockVectorSelector<const std::vector<VectorType>, false>
  {
    using BaseVectorType = const VectorType;
 
    static const BaseVectorType *
    get_vector_component(const std::vector<VectorType> &vec,
                         const unsigned int             component)
    {
      AssertIndexRange(component, vec.size());
      return &vec[component];
    }
  };
 
  template <typename VectorType>
  struct BlockVectorSelector<std::vector<VectorType *>, false>
  {
    using BaseVectorType = VectorType;
 
    static BaseVectorType *
    get_vector_component(std::vector<VectorType *> &vec,
                         const unsigned int         component)
    {
      AssertIndexRange(component, vec.size());
      return vec[component];
    }
  };
 
  template <typename VectorType>
  struct BlockVectorSelector<const std::vector<VectorType *>, false>
  {
    using BaseVectorType = const VectorType;
 
    static const BaseVectorType *
    get_vector_component(const std::vector<VectorType *> &vec,
                         const unsigned int               component)
    {
      AssertIndexRange(component, vec.size());
      return vec[component];
    }
  };
 
  template <typename VectorType, std::size_t N>
  struct BlockVectorSelector<std::array<VectorType *, N>, false>
  {
    using BaseVectorType = VectorType;
 
    static BaseVectorType *
    get_vector_component(std::array<VectorType *, N> &vec,
                         const unsigned int           component)
    {
      AssertIndexRange(component, vec.size());
      return vec[component];
    }
  };
} // namespace internal
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <typename VectorType, typename VectorOperation>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  read_write_operation(
    const VectorOperation                         &operation,
    const std::array<VectorType *, n_components_> &src,
    const std::array<
      const std::vector<ArrayView<const typename VectorType::value_type>> *,
      n_components_>           &src_sm,
    const std::bitset<n_lanes> &mask,
    const bool                  apply_constraints) const
{
  // Case 1: No MatrixFree object given, simple case because we do not need to
  // process constraints and need not care about vectorization -> go to
  // separate function
  if (this->matrix_free == nullptr)
    {
      read_write_operation_global(operation, src);
      return;
    }
 
  Assert(this->dof_info != nullptr, ExcNotInitialized());
  const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
  Assert(this->matrix_free->indices_initialized() == true, ExcNotInitialized());
  if (this->n_fe_components == 1)
    for (unsigned int comp = 0; comp < n_components; ++comp)
      {
        Assert(src[comp] != nullptr,
               ExcMessage("The finite element underlying this FEEvaluation "
                          "object is scalar, but you requested " +
                          std::to_string(n_components) +
                          " components via the template argument in "
                          "FEEvaluation. In that case, you must pass an "
                          "std::vector<VectorType> or a BlockVector to " +
                          "read_dof_values and distribute_local_to_global."));
        internal::check_vector_compatibility(*src[comp],
                                             *this->matrix_free,
                                             *this->dof_info);
      }
  else
    {
      internal::check_vector_compatibility(*src[0],
                                           *this->matrix_free,
                                           *this->dof_info);
    }
 
  const bool accesses_exterior_dofs =
    this->dof_access_index ==
      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
    this->is_interior_face() == false;
 
  // Case 2: contiguous indices which use reduced storage of indices and can
  // use vectorized load/store operations -> go to separate function
  if (this->cell != numbers::invalid_unsigned_int)
    {
      AssertIndexRange(
        this->cell,
        dof_info.index_storage_variants[this->dof_access_index].size());
 
      bool is_contiguous = true;
      // check if exterior cells are not contiguous (ECL case)
      if (accesses_exterior_dofs)
        {
          const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
          const unsigned int                       n_filled_lanes =
            dof_info.n_vectorization_lanes_filled
              [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
              [this->cell];
          // we have to check all filled lanes which are active in the mask
          for (unsigned int v = 0; v < n_filled_lanes; ++v)
            if (mask[v] == true &&
                dof_info.index_storage_variants
                    [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
                    [cells[v] / n_lanes] <
                  internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                    contiguous)
              is_contiguous = false;
        } // or if cell/face batch is not contiguous
      else if (dof_info.index_storage_variants
                 [is_face ?
                    this->dof_access_index :
                    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
                 [this->cell] < internal::MatrixFreeFunctions::DoFInfo::
                                  IndexStorageVariants::contiguous)
        {
          is_contiguous = false;
        }
 
      if (is_contiguous)
        {
          read_write_operation_contiguous(operation, src, src_sm, mask);
          return;
        }
    }
 
  // Case 3: standard operation with one index per degree of freedom -> go on
  // here
  std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
 
  const bool masking_is_active = mask.count() < n_lanes;
  if (masking_is_active)
    for (unsigned int v = 0; v < n_lanes; ++v)
      if (mask[v] == false)
        cells[v] = numbers::invalid_unsigned_int;
 
  std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
    vector_selector;
 
  const bool use_vectorized_path =
    !(masking_is_active || accesses_exterior_dofs);
 
  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
 
  if (apply_constraints && this->cell != numbers::invalid_unsigned_int &&
      dof_info.index_storage_variants
          [is_face ? this->dof_access_index :
                     internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
          [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
                            IndexStorageVariants::interleaved &&
      use_vectorized_path)
    {
      const unsigned int *dof_indices =
        dof_info.dof_indices_interleaved.data() +
        dof_info.row_starts[this->cell * this->n_fe_components * n_lanes]
          .first +
        this->dof_info
            ->component_dof_indices_offset[this->active_fe_index]
                                          [this->first_selected_component] *
          n_lanes;
 
      std::array<typename VectorType::value_type *, n_components> src_ptrs;
      if (n_components == 1 || this->n_fe_components == 1)
        for (unsigned int comp = 0; comp < n_components; ++comp)
          src_ptrs[comp] =
            const_cast<typename VectorType::value_type *>(src[comp]->begin());
      else
        src_ptrs[0] =
          const_cast<typename VectorType::value_type *>(src[0]->begin());
 
      if (n_components == 1 || this->n_fe_components == 1)
        for (unsigned int i = 0; i < dofs_per_component;
             ++i, dof_indices += n_lanes)
          for (unsigned int comp = 0; comp < n_components; ++comp)
            operation.process_dof_gather(
              dof_indices,
              *src[comp],
              0,
              src_ptrs[comp],
              this->values_dofs[comp * dofs_per_component + i],
              vector_selector);
      else
        for (unsigned int comp = 0; comp < n_components; ++comp)
          for (unsigned int i = 0; i < dofs_per_component;
               ++i, dof_indices += n_lanes)
            operation.process_dof_gather(
              dof_indices,
              *src[0],
              0,
              src_ptrs[0],
              this->values_dofs[comp * dofs_per_component + i],
              vector_selector);
      return;
    }
 
  std::array<VectorizedArrayType *, n_components> values_dofs;
  for (unsigned int c = 0; c < n_components; ++c)
    values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
                     c * dofs_per_component;
 
  // Allocate pointers, then initialize all of them to nullptrs and
  // below overwrite the ones we actually use:
  std::array<const unsigned int *, n_lanes> dof_indices;
  dof_indices.fill(nullptr);
 
  // Assign the appropriate cell ids for face/cell case and get the pointers
  // to the dof indices of the cells on all lanes
 
  bool               has_constraints = false;
  const unsigned int n_components_read =
    this->n_fe_components > 1 ? n_components : 1;
 
  if (is_face)
    {
      for (unsigned int v = 0; v < n_lanes; ++v)
        {
          if (cells[v] == numbers::invalid_unsigned_int)
            continue;
 
          Assert(cells[v] < dof_info.row_starts.size() - 1, ExcInternalError());
          const std::pair<unsigned int, unsigned int> *my_index_start =
            &dof_info.row_starts[cells[v] * this->n_fe_components +
                                 this->first_selected_component];
 
          // check whether any of the SIMD lanes has constraints, i.e., the
          // constraint indicator which is the second entry of row_starts
          // increments on this cell
          if (my_index_start[n_components_read].second !=
              my_index_start[0].second)
            has_constraints = true;
 
          dof_indices[v] =
            dof_info.dof_indices.data() + my_index_start[0].first;
        }
    }
  else
    {
      for (unsigned int v = 0; v < n_lanes; ++v)
        {
          if (cells[v] == numbers::invalid_unsigned_int)
            continue;
 
          const std::pair<unsigned int, unsigned int> *my_index_start =
            &dof_info.row_starts[cells[v] * this->n_fe_components +
                                 this->first_selected_component];
          if (my_index_start[n_components_read].second !=
              my_index_start[0].second)
            has_constraints = true;
 
          Assert(my_index_start[n_components_read].first ==
                     my_index_start[0].first ||
                   my_index_start[0].first < dof_info.dof_indices.size(),
                 ExcIndexRange(0,
                               my_index_start[0].first,
                               dof_info.dof_indices.size()));
          dof_indices[v] =
            dof_info.dof_indices.data() + my_index_start[0].first;
        }
    }
 
  for (unsigned int comp = 0; comp < n_components; ++comp)
    for (unsigned int i = 0; i < dofs_per_component; ++i)
      operation.process_empty(values_dofs[comp][i]);
 
  // Case where we have no constraints throughout the whole cell: Can go
  // through the list of DoFs directly
  if (!has_constraints && apply_constraints)
    {
      if (n_components == 1 || this->n_fe_components == 1)
        {
          for (unsigned int v = 0; v < n_lanes; ++v)
            {
              if (cells[v] == numbers::invalid_unsigned_int)
                continue;
 
              for (unsigned int i = 0; i < dofs_per_component; ++i)
                if (dof_indices[v][i] != numbers::invalid_unsigned_int)
                  for (unsigned int comp = 0; comp < n_components; ++comp)
                    operation.process_dof(dof_indices[v][i],
                                          *src[comp],
                                          values_dofs[comp][i][v]);
            }
        }
      else
        {
          for (unsigned int comp = 0; comp < n_components; ++comp)
            for (unsigned int v = 0; v < n_lanes; ++v)
              {
                if (cells[v] == numbers::invalid_unsigned_int)
                  continue;
 
                for (unsigned int i = 0; i < dofs_per_component; ++i)
                  if (dof_indices[v][comp * dofs_per_component + i] !=
                      numbers::invalid_unsigned_int)
                    operation.process_dof(
                      dof_indices[v][comp * dofs_per_component + i],
                      *src[0],
                      values_dofs[comp][i][v]);
              }
        }
      return;
    }
 
  // In the case where there are some constraints to be resolved, loop over
  // all vector components that are filled and then over local dofs. ind_local
  // holds local number on cell, index iterates over the elements of
  // index_local_to_global and dof_indices points to the global indices stored
  // in index_local_to_global
 
  for (unsigned int v = 0; v < n_lanes; ++v)
    {
      if (cells[v] == numbers::invalid_unsigned_int)
        continue;
 
      const unsigned int cell_index = cells[v];
      const unsigned int cell_dof_index =
        cell_index * this->n_fe_components + this->first_selected_component;
      unsigned int index_indicators =
        dof_info.row_starts[cell_dof_index].second;
      unsigned int next_index_indicators =
        dof_info.row_starts[cell_dof_index + 1].second;
 
      // For read_dof_values_plain, redirect the dof_indices field to the
      // unconstrained indices
      if (apply_constraints == false &&
          dof_info.row_starts_plain_indices[cell_index] !=
            numbers::invalid_unsigned_int)
        {
          dof_indices[v] =
            dof_info.plain_dof_indices.data() +
            this->dof_info
              ->component_dof_indices_offset[this->active_fe_index]
                                            [this->first_selected_component] +
            dof_info.row_starts_plain_indices[cell_index];
          next_index_indicators = index_indicators;
        }
 
      if (n_components == 1 || this->n_fe_components == 1)
        {
          unsigned int ind_local = 0;
          for (; index_indicators != next_index_indicators; ++index_indicators)
            {
              const std::pair<unsigned short, unsigned short> indicator =
                dof_info.constraint_indicator[index_indicators];
              // run through values up to next constraint
              for (unsigned int j = 0; j < indicator.first; ++j)
                if (dof_indices[v][j] != numbers::invalid_unsigned_int)
                  for (unsigned int comp = 0; comp < n_components; ++comp)
                    operation.process_dof(dof_indices[v][j],
                                          *src[comp],
                                          values_dofs[comp][ind_local + j][v]);
 
              ind_local += indicator.first;
              dof_indices[v] += indicator.first;
 
              // constrained case: build the local value as a linear
              // combination of the global value according to constraints
              Number value[n_components];
              for (unsigned int comp = 0; comp < n_components; ++comp)
                operation.pre_constraints(values_dofs[comp][ind_local][v],
                                          value[comp]);
 
              const Number *data_val =
                this->matrix_free->constraint_pool_begin(indicator.second);
              const Number *end_pool =
                this->matrix_free->constraint_pool_end(indicator.second);
              for (; data_val != end_pool; ++data_val, ++dof_indices[v])
                if (*dof_indices[v] != numbers::invalid_unsigned_int)
                  for (unsigned int comp = 0; comp < n_components; ++comp)
                    operation.process_constraint(*dof_indices[v],
                                                 *data_val,
                                                 *src[comp],
                                                 value[comp]);
 
              for (unsigned int comp = 0; comp < n_components; ++comp)
                operation.post_constraints(value[comp],
                                           values_dofs[comp][ind_local][v]);
              ++ind_local;
            }
 
          AssertIndexRange(ind_local, dofs_per_component + 1);
 
          for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
            if (*dof_indices[v] != numbers::invalid_unsigned_int)
              for (unsigned int comp = 0; comp < n_components; ++comp)
                operation.process_dof(*dof_indices[v],
                                      *src[comp],
                                      values_dofs[comp][ind_local][v]);
        }
      else
        {
          // case with vector-valued finite elements where all components are
          // included in one single vector. Assumption: first come all entries
          // to the first component, then all entries to the second one, and
          // so on. This is ensured by the way MatrixFree reads out the
          // indices.
          for (unsigned int comp = 0; comp < n_components; ++comp)
            {
              unsigned int ind_local = 0;
 
              // check whether there is any constraint on the current cell
              for (; index_indicators != next_index_indicators;
                   ++index_indicators)
                {
                  const std::pair<unsigned short, unsigned short> indicator =
                    dof_info.constraint_indicator[index_indicators];
 
                  // run through values up to next constraint
                  for (unsigned int j = 0; j < indicator.first; ++j)
                    if (dof_indices[v][j] != numbers::invalid_unsigned_int)
                      operation.process_dof(
                        dof_indices[v][j],
                        *src[0],
                        values_dofs[comp][ind_local + j][v]);
                  ind_local += indicator.first;
                  dof_indices[v] += indicator.first;
 
                  // constrained case: build the local value as a linear
                  // combination of the global value according to constraints
                  Number value;
                  operation.pre_constraints(values_dofs[comp][ind_local][v],
                                            value);
 
                  const Number *data_val =
                    this->matrix_free->constraint_pool_begin(indicator.second);
                  const Number *end_pool =
                    this->matrix_free->constraint_pool_end(indicator.second);
 
                  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
                    if (*dof_indices[v] != numbers::invalid_unsigned_int)
                      operation.process_constraint(*dof_indices[v],
                                                   *data_val,
                                                   *src[0],
                                                   value);
 
                  operation.post_constraints(value,
                                             values_dofs[comp][ind_local][v]);
                  ++ind_local;
                }
 
              AssertIndexRange(ind_local, dofs_per_component + 1);
 
              // get the dof values past the last constraint
              for (; ind_local < dofs_per_component;
                   ++dof_indices[v], ++ind_local)
                if (*dof_indices[v] != numbers::invalid_unsigned_int)
                  {
                    AssertIndexRange(*dof_indices[v], src[0]->size());
                    operation.process_dof(*dof_indices[v],
                                          *src[0],
                                          values_dofs[comp][ind_local][v]);
                  }
 
              if (apply_constraints == true && comp + 1 < n_components)
                next_index_indicators =
                  dof_info.row_starts[cell_dof_index + comp + 2].second;
            }
        }
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <typename VectorType, typename VectorOperation>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  read_write_operation_global(
    const VectorOperation                         &operation,
    const std::array<VectorType *, n_components_> &src) const
{
  Assert(!local_dof_indices.empty(), ExcNotInitialized());
 
  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
  unsigned int      index = this->first_selected_component * dofs_per_component;
  for (unsigned int comp = 0; comp < n_components; ++comp)
    {
      for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
        {
          operation.process_empty(
            this->values_dofs[comp * dofs_per_component + i]);
          operation.process_dof_global(
            local_dof_indices[this->data->lexicographic_numbering[index]],
            *src[0],
            this->values_dofs[comp * dofs_per_component + i][0]);
        }
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <typename VectorType, typename VectorOperation>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  read_write_operation_contiguous(
    const VectorOperation                         &operation,
    const std::array<VectorType *, n_components_> &src,
    const std::array<
      const std::vector<ArrayView<const typename VectorType::value_type>> *,
      n_components_>           &vectors_sm,
    const std::bitset<n_lanes> &mask) const
{
  // This functions processes the functions read_dof_values,
  // distribute_local_to_global, and set_dof_values with the same code for
  // contiguous cell indices (DG case). The distinction between these three
  // cases is made by the input VectorOperation that either reads values from
  // a vector and puts the data into the local data field or write local data
  // into the vector. Certain operations are no-ops for the given use case.
 
  std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
                                                               vector_selector;
  const internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex ind =
    is_face ? this->dof_access_index :
              internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
  const unsigned int n_active_lanes = mask.count();
 
  const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
  const std::vector<unsigned int>              &dof_indices_cont =
    dof_info.dof_indices_contiguous[ind];
 
  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
  std::array<VectorizedArrayType *, n_components> values_dofs{{nullptr}};
  for (unsigned int c = 0; c < n_components; ++c)
    values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
                     c * dofs_per_component;
 
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotImplemented());
 
  const bool accesses_exterior_dofs =
    this->dof_access_index ==
      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
    this->is_interior_face() == false;
 
  // Whether the n_components copies of a vector-valued quantity live in
  // n_components separate (typically scalar) vectors, one component per
  // vector -- as opposed to all components being interleaved within a
  // single vector that belongs to an FESystem with n_fe_components ==
  // n_components. This distinction recurs throughout this function
  // because the two cases require different index computations; we
  // evaluate it once here rather than at each use.
  const bool separate_vectors =
    (n_components == 1 || this->n_fe_components == 1);
 
  // Simple case: We have contiguous storage, so we can simply copy out the
  // data
  if (dof_info.index_storage_variants[ind][this->cell] ==
        internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
          interleaved_contiguous &&
      n_active_lanes == n_lanes && !accesses_exterior_dofs)
    {
      const unsigned int dof_index =
        dof_indices_cont[this->cell * n_lanes] +
        this->dof_info
            ->component_dof_indices_offset[this->active_fe_index]
                                          [this->first_selected_component] *
          n_lanes;
      if (separate_vectors)
        for (unsigned int comp = 0; comp < n_components; ++comp)
          operation.process_dofs_vectorized(dofs_per_component,
                                            dof_index,
                                            *src[comp],
                                            values_dofs[comp],
                                            vector_selector);
      else
        operation.process_dofs_vectorized(dofs_per_component * n_components,
                                          dof_index,
                                          *src[0],
                                          values_dofs[0],
                                          vector_selector);
      return;
    }
 
  const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
 
  // More general case: Must go through the components one by one and apply
  // some transformations
  const unsigned int n_filled_lanes =
    dof_info.n_vectorization_lanes_filled[ind][this->cell];
 
  const bool use_vectorized_path = n_filled_lanes == n_lanes &&
                                   n_active_lanes == n_lanes &&
                                   !accesses_exterior_dofs;
 
  if (vectors_sm[0] != nullptr)
    {
      const auto compute_vector_ptrs = [&](const unsigned int comp) {
        std::array<typename VectorType::value_type *, n_lanes> vector_ptrs{
          {nullptr}};
 
        const auto upper_bound =
          std::min<unsigned int>(n_filled_lanes, n_lanes);
        for (unsigned int v = 0; v < upper_bound; ++v)
          {
            if (mask[v] == false)
              {
                vector_ptrs[v] = nullptr;
                continue;
              }
 
            Assert(cells[v] != numbers::invalid_unsigned_int,
                   ExcNotImplemented());
            Assert(ind < dof_info.dof_indices_contiguous_sm.size(),
                   ExcIndexRange(ind,
                                 0,
                                 dof_info.dof_indices_contiguous_sm.size()));
            Assert(
              cells[v] < dof_info.dof_indices_contiguous_sm[ind].size(),
              ExcIndexRange(cells[v],
                            0,
                            dof_info.dof_indices_contiguous_sm[ind].size()));
 
            const auto &temp =
              dof_info.dof_indices_contiguous_sm[ind][cells[v]];
 
            if (temp.first != numbers::invalid_unsigned_int)
              vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
                vectors_sm[comp]->operator[](temp.first).data() + temp.second +
                dof_info.component_dof_indices_offset
                  [this->active_fe_index][this->first_selected_component]);
            else
              vector_ptrs[v] = nullptr;
          }
        for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
          vector_ptrs[v] = nullptr;
 
        return vector_ptrs;
      };
 
      if (use_vectorized_path)
        {
          // If the components live in separate vectors, gather each one
          // with its own transposed access. If instead all components
          // come from a single FESystem-based vector, their DoFs are
          // contiguous per cell in both the source vector and
          // values_dofs, so we can process all of them with a single
          // bulk transposed access instead -- one loop iteration with a
          // correspondingly larger block size.
          const unsigned int n_blocks   = separate_vectors ? n_components : 1;
          const unsigned int block_size = separate_vectors ?
                                            dofs_per_component :
                                            dofs_per_component * n_components;
          for (unsigned int comp = 0; comp < n_blocks; ++comp)
            {
              auto vector_ptrs =
                compute_vector_ptrs(separate_vectors ? comp : 0);
              operation.process_dofs_vectorized_transpose(block_size,
                                                          vector_ptrs,
                                                          separate_vectors ?
                                                            values_dofs[comp] :
                                                            &values_dofs[0][0],
                                                          vector_selector);
            }
        }
      else
        for (unsigned int comp = 0; comp < n_components; ++comp)
          {
            auto vector_ptrs = compute_vector_ptrs(separate_vectors ? comp : 0);
 
            for (unsigned int i = 0; i < dofs_per_component; ++i)
              operation.process_empty(values_dofs[comp][i]);
 
            // As above: with separate vectors, component comp's DoFs sit
            // at the start of vector_ptrs; with a single FESystem-based
            // vector, they are offset by the preceding components' DoFs.
            const unsigned int comp_offset =
              separate_vectors ? 0 : comp * dofs_per_component;
            for (unsigned int v = 0; v < n_filled_lanes; ++v)
              if (mask[v] == true)
                for (unsigned int i = 0; i < dofs_per_component; ++i)
                  operation.process_dof(vector_ptrs[v][i + comp_offset],
                                        values_dofs[comp][i][v]);
          }
      return;
    }
 
  std::array<unsigned int, n_lanes> dof_indices{
    {numbers::invalid_unsigned_int}};
  Assert(n_filled_lanes <= n_lanes, ExcInternalError());
  for (unsigned int v = 0; v < n_filled_lanes; ++v)
    {
      Assert(mask[v] == false || cells[v] != numbers::invalid_unsigned_int,
             ExcNotImplemented());
      if (mask[v] == true)
        dof_indices[v] =
          dof_indices_cont[cells[v]] +
          this->dof_info
              ->component_dof_indices_offset[this->active_fe_index]
                                            [this->first_selected_component] *
            dof_info.dof_indices_interleave_strides[ind][cells[v]];
    }
 
  // In the case with contiguous cell indices, we know that there are no
  // constraints and that the indices within each element are contiguous
  if (use_vectorized_path)
    {
      if (dof_info.index_storage_variants[ind][this->cell] ==
          internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
            contiguous)
        {
          if (separate_vectors)
            for (unsigned int comp = 0; comp < n_components; ++comp)
              operation.process_dofs_vectorized_transpose(dofs_per_component,
                                                          dof_indices.data(),
                                                          *src[comp],
                                                          values_dofs[comp],
                                                          vector_selector);
          else
            operation.process_dofs_vectorized_transpose(dofs_per_component *
                                                          n_components,
                                                        dof_indices.data(),
                                                        *src[0],
                                                        &values_dofs[0][0],
                                                        vector_selector);
        }
      else if (dof_info.index_storage_variants[ind][this->cell] ==
               internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                 interleaved_contiguous_strided)
        {
          std::array<typename VectorType::value_type *, n_components> src_ptrs{
            {nullptr}};
          if (separate_vectors)
            for (unsigned int comp = 0; comp < n_components; ++comp)
              src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
                src[comp]->begin());
          else
            src_ptrs[0] =
              const_cast<typename VectorType::value_type *>(src[0]->begin());
 
          if (separate_vectors)
            for (unsigned int i = 0; i < dofs_per_component; ++i)
              {
                for (unsigned int comp = 0; comp < n_components; ++comp)
                  operation.process_dof_gather(dof_indices.data(),
                                               *src[comp],
                                               i * n_lanes,
                                               src_ptrs[comp] + i * n_lanes,
                                               values_dofs[comp][i],
                                               vector_selector);
              }
          else
            for (unsigned int comp = 0; comp < n_components; ++comp)
              for (unsigned int i = 0; i < dofs_per_component; ++i)
                {
                  operation.process_dof_gather(
                    dof_indices.data(),
                    *src[0],
                    (comp * dofs_per_component + i) * n_lanes,
                    src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
                    values_dofs[comp][i],
                    vector_selector);
                }
        }
      else
        {
          Assert(dof_info.index_storage_variants[ind][this->cell] ==
                   internal::MatrixFreeFunctions::DoFInfo::
                     IndexStorageVariants::interleaved_contiguous_mixed_strides,
                 ExcNotImplemented());
          std::array<typename VectorType::value_type *, n_components> src_ptrs{
            {nullptr}};
          if (separate_vectors)
            for (unsigned int comp = 0; comp < n_components; ++comp)
              src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
                src[comp]->begin());
          else
            src_ptrs[0] =
              const_cast<typename VectorType::value_type *>(src[0]->begin());
 
          const unsigned int *offsets =
            &dof_info.dof_indices_interleave_strides[ind][n_lanes * this->cell];
          if (separate_vectors)
            for (unsigned int i = 0; i < dofs_per_component; ++i)
              {
                for (unsigned int comp = 0; comp < n_components; ++comp)
                  operation.process_dof_gather(dof_indices.data(),
                                               *src[comp],
                                               0,
                                               src_ptrs[comp],
                                               values_dofs[comp][i],
                                               vector_selector);
                DEAL_II_OPENMP_SIMD_PRAGMA
                for (unsigned int v = 0; v < n_lanes; ++v)
                  dof_indices[v] += offsets[v];
              }
          else
            for (unsigned int comp = 0; comp < n_components; ++comp)
              for (unsigned int i = 0; i < dofs_per_component; ++i)
                {
                  operation.process_dof_gather(dof_indices.data(),
                                               *src[0],
                                               0,
                                               src_ptrs[0],
                                               values_dofs[comp][i],
                                               vector_selector);
                  DEAL_II_OPENMP_SIMD_PRAGMA
                  for (unsigned int v = 0; v < n_lanes; ++v)
                    dof_indices[v] += offsets[v];
                }
        }
    }
  else
    for (unsigned int comp = 0; comp < n_components; ++comp)
      {
        for (unsigned int i = 0; i < dofs_per_component; ++i)
          operation.process_empty(values_dofs[comp][i]);
        if (accesses_exterior_dofs)
          {
            for (unsigned int v = 0; v < n_filled_lanes; ++v)
              if (mask[v] == true)
                {
                  if (dof_info.index_storage_variants
                        [ind][cells[v] / VectorizedArrayType::size()] ==
                      internal::MatrixFreeFunctions::DoFInfo::
                        IndexStorageVariants::contiguous)
                    {
                      if (separate_vectors)
                        {
                          for (unsigned int i = 0; i < dofs_per_component; ++i)
                            operation.process_dof(dof_indices[v] + i,
                                                  *src[comp],
                                                  values_dofs[comp][i][v]);
                        }
                      else
                        {
                          for (unsigned int i = 0; i < dofs_per_component; ++i)
                            operation.process_dof(dof_indices[v] + i +
                                                    comp * dofs_per_component,
                                                  *src[0],
                                                  values_dofs[comp][i][v]);
                        }
                    }
                  else
                    {
                      const unsigned int offset =
                        dof_info.dof_indices_interleave_strides[ind][cells[v]];
                      AssertIndexRange(offset, VectorizedArrayType::size() + 1);
                      if (separate_vectors)
                        {
                          for (unsigned int i = 0; i < dofs_per_component; ++i)
                            operation.process_dof(dof_indices[v] + i * offset,
                                                  *src[comp],
                                                  values_dofs[comp][i][v]);
                        }
                      else
                        {
                          for (unsigned int i = 0; i < dofs_per_component; ++i)
                            operation.process_dof(
                              dof_indices[v] +
                                (i + comp * dofs_per_component) * offset,
                              *src[0],
                              values_dofs[comp][i][v]);
                        }
                    }
                }
          }
        else
          {
            if (dof_info.index_storage_variants[ind][this->cell] ==
                internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                  contiguous)
              {
                if (separate_vectors)
                  {
                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
                      if (mask[v] == true)
                        for (unsigned int i = 0; i < dofs_per_component; ++i)
                          operation.process_dof(dof_indices[v] + i,
                                                *src[comp],
                                                values_dofs[comp][i][v]);
                  }
                else
                  {
                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
                      if (mask[v] == true)
                        for (unsigned int i = 0; i < dofs_per_component; ++i)
                          operation.process_dof(dof_indices[v] + i +
                                                  comp * dofs_per_component,
                                                *src[0],
                                                values_dofs[comp][i][v]);
                  }
              }
            else
              {
                const unsigned int *offsets =
                  &dof_info.dof_indices_interleave_strides
                     [ind][VectorizedArrayType::size() * this->cell];
                for (unsigned int v = 0; v < n_filled_lanes; ++v)
                  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
                if (separate_vectors)
                  for (unsigned int v = 0; v < n_filled_lanes; ++v)
                    {
                      if (mask[v] == true)
                        for (unsigned int i = 0; i < dofs_per_component; ++i)
                          operation.process_dof(dof_indices[v] + i * offsets[v],
                                                *src[comp],
                                                values_dofs[comp][i][v]);
                    }
                else
                  {
                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
                      if (mask[v] == true)
                        for (unsigned int i = 0; i < dofs_per_component; ++i)
                          operation.process_dof(
                            dof_indices[v] +
                              (i + comp * dofs_per_component) * offsets[v],
                            *src[0],
                            values_dofs[comp][i][v]);
                  }
              }
          }
      }
}
 
namespace internal
{
  template <
    typename Number,
    typename VectorType,
    std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * = nullptr>
  decltype(std::declval<VectorType>().begin())
  get_beginning(VectorType &vec)
  {
    return vec.begin();
  }
 
  template <
    typename Number,
    typename VectorType,
    std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * = nullptr>
  typename VectorType::value_type *
  get_beginning(VectorType &)
  {
    return nullptr;
  }
 
  template <typename VectorType,
            std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
              nullptr>
  const std::vector<ArrayView<const typename VectorType::value_type>> *
  get_shared_vector_data(VectorType        *vec,
                         const bool         is_valid_mode_for_sm,
                         const unsigned int active_fe_index,
                         const internal::MatrixFreeFunctions::DoFInfo *dof_info)
  {
    // note: no hp is supported
    if (is_valid_mode_for_sm &&
        dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
            .size() > 0 &&
        active_fe_index == 0)
      return &vec->shared_vector_data();
    else
      return nullptr;
  }
 
  template <typename VectorType,
            std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
              * = nullptr>
  const std::vector<ArrayView<const typename VectorType::value_type>> *
  get_shared_vector_data(VectorType *,
                         const bool,
                         const unsigned int,
                         const internal::MatrixFreeFunctions::DoFInfo *)
  {
    return nullptr;
  }
 
  template <int n_components, typename VectorType>
  std::pair<
    std::array<typename internal::BlockVectorSelector<
                 VectorType,
                 IsBlockVector<VectorType>::value>::BaseVectorType *,
               n_components>,
    std::array<
      const std::vector<ArrayView<const typename internal::BlockVectorSelector<
        VectorType,
        IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
      n_components>>
  get_vector_data(VectorType        &src,
                  const unsigned int first_index,
                  const bool         is_valid_mode_for_sm,
                  const unsigned int active_fe_index,
                  const internal::MatrixFreeFunctions::DoFInfo *dof_info)
  {
    // select between block vectors and non-block vectors. Note that the number
    // of components is checked in the internal data
    std::pair<
      std::array<typename internal::BlockVectorSelector<
                   VectorType,
                   IsBlockVector<VectorType>::value>::BaseVectorType *,
                 n_components>,
      std::array<
        const std::vector<
          ArrayView<const typename internal::BlockVectorSelector<
            VectorType,
            IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
        n_components>>
      src_data;
 
    for (unsigned int d = 0; d < n_components; ++d)
      src_data.first[d] = internal::BlockVectorSelector<
        VectorType,
        IsBlockVector<VectorType>::value>::get_vector_component(src,
                                                                d +
                                                                  first_index);
 
    for (unsigned int d = 0; d < n_components; ++d)
      src_data.second[d] = get_shared_vector_data(
        const_cast<typename internal::BlockVectorSelector<
          std::remove_const_t<VectorType>,
          IsBlockVector<std::remove_const_t<VectorType>>::value>::BaseVectorType
                     *>(src_data.first[d]),
        is_valid_mode_for_sm,
        active_fe_index,
        dof_info);
 
    return src_data;
  }
} // namespace internal
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  apply_hanging_node_constraints(const bool transpose) const
{
  if (this->dof_info == nullptr ||
      this->dof_info->hanging_node_constraint_masks.empty() ||
      this->dof_info->hanging_node_constraint_masks_comp.empty() ||
      this->dof_info->hanging_node_constraint_masks_comp
          [this->active_fe_index][this->first_selected_component] == false)
    return; // nothing to do with faces
 
  std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
    constraint_mask{{internal::MatrixFreeFunctions::
                       unconstrained_compressed_constraint_kind}};
 
  bool hn_available = false;
 
  const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
 
  for (unsigned int v = 0; v < n_lanes; ++v)
    {
      if (cells[v] == numbers::invalid_unsigned_int)
        {
          constraint_mask[v] = internal::MatrixFreeFunctions::
            unconstrained_compressed_constraint_kind;
          continue;
        }
 
      const unsigned int cell_index = cells[v];
      const auto         mask =
        this->dof_info->hanging_node_constraint_masks[cell_index];
      constraint_mask[v] = mask;
 
      hn_available |= (mask != internal::MatrixFreeFunctions::
                                 unconstrained_compressed_constraint_kind);
    }
 
  if (hn_available == false)
    return; // no hanging node on cell batch -> nothing to do
 
  internal::FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::
    apply(n_components,
          this->data->data.front().fe_degree,
          this->get_shape_info(),
          transpose,
          constraint_mask,
          this->values_dofs);
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  read_dof_values(const VectorType           &src,
                  const unsigned int          first_index,
                  const std::bitset<n_lanes> &mask)
{
  const auto src_data = internal::get_vector_data<n_components_>(
    src,
    first_index,
    this->dof_info != nullptr &&
      this->dof_access_index ==
        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
    this->active_fe_index,
    this->dof_info);
 
  internal::VectorReader<Number, VectorizedArrayType> reader;
  read_write_operation(reader, src_data.first, src_data.second, mask, true);
 
  apply_hanging_node_constraints(false);
 
  if constexpr (running_in_debug_mode())
    {
      this->dof_values_initialized = true;
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  read_dof_values_plain(const VectorType           &src,
                        const unsigned int          first_index,
                        const std::bitset<n_lanes> &mask)
{
  const auto src_data = internal::get_vector_data<n_components_>(
    src,
    first_index,
    this->dof_info != nullptr &&
      this->dof_access_index ==
        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
    this->active_fe_index,
    this->dof_info);
 
  internal::VectorReader<Number, VectorizedArrayType> reader;
  read_write_operation(reader, src_data.first, src_data.second, mask, false);
 
  if constexpr (running_in_debug_mode())
    {
      this->dof_values_initialized = true;
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  distribute_local_to_global(VectorType                 &dst,
                             const unsigned int          first_index,
                             const std::bitset<n_lanes> &mask) const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->dof_values_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
 
  apply_hanging_node_constraints(true);
 
  const auto dst_data = internal::get_vector_data<n_components_>(
    dst,
    first_index,
    this->dof_info != nullptr &&
      this->dof_access_index ==
        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
    this->active_fe_index,
    this->dof_info);
 
  internal::VectorDistributorLocalToGlobal<Number, VectorizedArrayType>
    distributor;
  read_write_operation(distributor, dst_data.first, dst_data.second, mask);
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  set_dof_values(VectorType                 &dst,
                 const unsigned int          first_index,
                 const std::bitset<n_lanes> &mask) const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->dof_values_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
 
  const auto dst_data = internal::get_vector_data<n_components_>(
    dst,
    first_index,
    this->dof_info != nullptr &&
      this->dof_access_index ==
        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
    this->active_fe_index,
    this->dof_info);
 
  internal::VectorSetter<Number, VectorizedArrayType> setter;
  read_write_operation(setter, dst_data.first, dst_data.second, mask);
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  set_dof_values_plain(VectorType                 &dst,
                       const unsigned int          first_index,
                       const std::bitset<n_lanes> &mask) const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->dof_values_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
 
  const auto dst_data = internal::get_vector_data<n_components_>(
    dst,
    first_index,
    this->dof_access_index ==
      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
    this->active_fe_index,
    this->dof_info);
 
  internal::VectorSetter<Number, VectorizedArrayType> setter;
  read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
}
 
 
 
/*------------------------------ access to data fields ----------------------*/
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE
  typename FEEvaluationBase<dim,
                            n_components_,
                            Number,
                            is_face,
                            VectorizedArrayType>::value_type
  FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
    get_dof_value(const unsigned int dof) const
{
  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
  if constexpr (n_components == 1)
    return this->values_dofs[dof];
  else
    {
      const std::size_t dofs = this->data->dofs_per_component_on_cell;
      Tensor<1, n_components_, VectorizedArrayType> return_value;
      for (unsigned int comp = 0; comp < n_components; ++comp)
        return_value[comp] = this->values_dofs[comp * dofs + dof];
      return return_value;
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE
  typename FEEvaluationBase<dim,
                            n_components_,
                            Number,
                            is_face,
                            VectorizedArrayType>::value_type
  FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
    get_value(const unsigned int q_point) const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->values_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
 
  AssertIndexRange(q_point, this->n_quadrature_points);
  if constexpr (n_components == 1)
    return this->values_quad[q_point];
  else
    {
      if (n_components == dim &&
          this->data->element_type ==
            internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
        {
          // Piola transform is required
          if constexpr (running_in_debug_mode())
            {
              Assert(this->values_quad_initialized == true,
                     internal::ExcAccessToUninitializedField());
            }
 
          AssertIndexRange(q_point, this->n_quadrature_points);
          Assert(this->J_value != nullptr,
                 internal::ExcMatrixFreeAccessToUninitializedMappingField(
                   "update_values"));
          const std::size_t nqp = this->n_quadrature_points;
          Tensor<1, n_components, VectorizedArrayType> value_out;
 
          if (!is_face &&
              this->cell_type == internal::MatrixFreeFunctions::cartesian)
            {
              // Cartesian cell
              const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
              const VectorizedArrayType                 inv_det =
                (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
                                             this->jacobian[0][0][0] * this->jacobian[0][1][1] *
                               this->jacobian[0][2][2];
 
              // J * u * det(J^-1)
              for (unsigned int comp = 0; comp < n_components; ++comp)
                value_out[comp] = this->values_quad[comp * nqp + q_point] *
                                  jac[comp][comp] * inv_det;
            }
          else
            {
              // Affine or general cell
              const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
                (this->cell_type > internal::MatrixFreeFunctions::affine) ?
                  this->jacobian[q_point] :
                  this->jacobian[0];
              const Tensor<2, dim, VectorizedArrayType> jac =
                (this->cell_type > internal::MatrixFreeFunctions::affine) ?
                  transpose(invert(inv_t_jac)) :
                  this->jacobian[1];
 
              // Derivatives are reordered for faces. Need to take this into
              // account
              const VectorizedArrayType inv_det =
                (is_face && dim == 2 && this->get_face_no() < 2) ?
                  -determinant(inv_t_jac) :
                  determinant(inv_t_jac);
              // J * u * det(J^-1)
              for (unsigned int comp = 0; comp < n_components; ++comp)
                {
                  value_out[comp] = this->values_quad[q_point] * jac[comp][0];
                  for (unsigned int e = 1; e < dim; ++e)
                    value_out[comp] +=
                      this->values_quad[e * nqp + q_point] * jac[comp][e];
                  value_out[comp] *= inv_det;
                }
            }
          return value_out;
        }
      else if (n_components == dim &&
               this->data->element_type ==
                 internal::MatrixFreeFunctions::ElementType::tensor_nedelec)
        {
          // Piola transform is required
          if constexpr (running_in_debug_mode())
            {
              Assert(this->values_quad_initialized == true,
                     internal::ExcAccessToUninitializedField());
            }
 
          AssertIndexRange(q_point, this->n_quadrature_points);
          Assert(this->J_value != nullptr,
                 internal::ExcMatrixFreeAccessToUninitializedMappingField(
                   "update_values"));
          const std::size_t nqp = this->n_quadrature_points;
          Tensor<1, n_components, VectorizedArrayType> value_out;
 
          if (!is_face &&
              this->cell_type == internal::MatrixFreeFunctions::cartesian)
            {
              // Cartesian cell
              const Tensor<2, dim, VectorizedArrayType> inv_transp_jac =
                this->jacobian[0];
 
 
              // J^{-T} * u
              for (unsigned int comp = 0; comp < n_components; ++comp)
                value_out[comp] = this->values_quad[comp * nqp + q_point] *
                                  inv_transp_jac[comp][comp];
            }
          else
            {
              // Affine or general cell
              const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
                (this->cell_type > internal::MatrixFreeFunctions::affine) ?
                  this->jacobian[q_point] :
                  this->jacobian[0];
 
              // J^{-T} * u
              for (unsigned int comp = 0; comp < n_components; ++comp)
                {
                  value_out[comp] =
                    this->values_quad[q_point] * inv_t_jac[comp][0];
                  for (unsigned int e = 1; e < dim; ++e)
                    value_out[comp] +=
                      this->values_quad[e * nqp + q_point] * inv_t_jac[comp][e];
                }
            }
          return value_out;
        }
      else
        {
          const std::size_t nqp = this->n_quadrature_points;
          Tensor<1, n_components_, VectorizedArrayType> return_value;
          for (unsigned int comp = 0; comp < n_components; ++comp)
            return_value[comp] = this->values_quad[comp * nqp + q_point];
          return return_value;
        }
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE
  typename FEEvaluationBase<dim,
                            n_components_,
                            Number,
                            is_face,
                            VectorizedArrayType>::gradient_type
  FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
    get_gradient(const unsigned int q_point) const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->gradients_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
 
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  const std::size_t nqp = this->n_quadrature_points;
 
  if constexpr (n_components == dim && dim > 1)
    {
      if (this->data->element_type ==
          internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
        {
          // Piola transform is required
          if constexpr (running_in_debug_mode())
            {
              Assert(this->gradients_quad_initialized == true,
                     internal::ExcAccessToUninitializedField());
            }
 
          AssertIndexRange(q_point, this->n_quadrature_points);
          Assert(this->jacobian != nullptr,
                 internal::ExcMatrixFreeAccessToUninitializedMappingField(
                   "update_gradients"));
          const std::size_t nqp   = this->n_quadrature_points;
          const std::size_t nqp_d = nqp * dim;
          Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_out;
          const VectorizedArrayType                          *gradients =
            this->gradients_quad + q_point * dim;
 
 
          if (!is_face &&
              this->cell_type == internal::MatrixFreeFunctions::cartesian)
            {
              // Cartesian cell
              const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
                this->jacobian[0];
              const Tensor<2, dim, VectorizedArrayType> &jac =
                this->jacobian[1];
              const VectorizedArrayType inv_det =
                (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
                             this->jacobian[0][0][0] * this->jacobian[0][1][1] *
                               this->jacobian[0][2][2];
 
              // J * grad_quad * J^-1 * det(J^-1)
              for (unsigned int d = 0; d < dim; ++d)
                for (unsigned int comp = 0; comp < n_components; ++comp)
                  grad_out[comp][d] = gradients[comp * nqp_d + d] *
                                      inv_t_jac[d][d] *
                                      (jac[comp][comp] * inv_det);
            }
          else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
            {
              // Affine cell
              const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
                this->jacobian[0];
              const Tensor<2, dim, VectorizedArrayType> &jac =
                this->jacobian[1];
 
              // Derivatives are reordered for faces. Need to take this into
              // account
              const VectorizedArrayType inv_det =
                (is_face && dim == 2 && this->get_face_no() < 2) ?
                  -determinant(inv_t_jac) :
                  determinant(inv_t_jac);
 
              VectorizedArrayType tmp[dim][dim];
              // J * grad_quad * J^-1 * det(J^-1)
              for (unsigned int d = 0; d < dim; ++d)
                for (unsigned int e = 0; e < dim; ++e)
                  {
                    tmp[d][e] = inv_t_jac[d][0] * gradients[e * nqp_d + 0];
                    for (unsigned int f = 1; f < dim; ++f)
                      tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
                  }
              for (unsigned int comp = 0; comp < n_components; ++comp)
                for (unsigned int d = 0; d < dim; ++d)
                  {
                    VectorizedArrayType res = jac[comp][0] * tmp[d][0];
                    for (unsigned int f = 1; f < dim; ++f)
                      res += jac[comp][f] * tmp[d][f];
 
                    grad_out[comp][d] = res * inv_det;
                  }
            }
          else
            {
              // General cell
 
              // This assert could be removed if we make sure that this is
              // updated even though update_hessians or update_jacobian_grads is
              // not passed, i.e make the necessary changes in
              // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags
              Assert(this->jacobian_gradients_non_inverse != nullptr,
                     internal::ExcMatrixFreeAccessToUninitializedMappingField(
                       "update_hessians"));
 
              const auto jac_grad =
                this->jacobian_gradients_non_inverse[q_point];
              const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
                this->jacobian[q_point];
 
              // Derivatives are reordered for faces. Need to take this into
              // account
              const VectorizedArrayType inv_det =
                (is_face && dim == 2 && this->get_face_no() < 2) ?
                  -determinant(inv_t_jac) :
                  determinant(inv_t_jac);
              const Tensor<2, dim, VectorizedArrayType> t_jac =
                invert(inv_t_jac);
 
              // (J * grad_quad) * J^-1 * det(J^-1), part in braces
              VectorizedArrayType tmp[dim][dim];
              for (unsigned int d = 0; d < dim; ++d)
                for (unsigned int e = 0; e < dim; ++e)
                  {
                    tmp[e][d] = t_jac[0][d] * gradients[0 * nqp_d + e];
                    for (unsigned int f = 1; f < dim; ++f)
                      tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
                  }
 
              // Add (jac_grad * values) * J^{-1} * det(J^{-1}), combine terms
              // outside braces with gradient part from above
              for (unsigned int d = 0; d < dim; ++d)
                {
                  for (unsigned int e = 0; e < dim; ++e)
                    tmp[e][d] +=
                      jac_grad[e][d] * this->values_quad[e * nqp + q_point];
                  for (unsigned int f = 0, r = dim; f < dim; ++f)
                    for (unsigned int k = f + 1; k < dim; ++k, ++r)
                      {
                        tmp[k][d] +=
                          jac_grad[r][d] * this->values_quad[f * nqp + q_point];
                        tmp[f][d] +=
                          jac_grad[r][d] * this->values_quad[k * nqp + q_point];
                      }
                }
 
              // Apply J^{-1} appearing in both terms outside braces above
              for (unsigned int d = 0; d < dim; ++d)
                for (unsigned int e = 0; e < dim; ++e)
                  {
                    VectorizedArrayType res = tmp[0][d] * inv_t_jac[e][0];
                    for (unsigned int f = 1; f < dim; ++f)
                      res += tmp[f][d] * inv_t_jac[e][f];
                    grad_out[d][e] = res;
                  }
 
              // Add -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1})),
              // which can be expressed as a rank-1 update tmp[d] * tmp4[e],
              // where tmp = J * values and tmp4 = (J^{-T} * jac_grad * J^{-1})
              VectorizedArrayType tmp3[dim], tmp4[dim];
              for (unsigned int d = 0; d < dim; ++d)
                {
                  tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
                  for (unsigned int e = 1; e < dim; ++e)
                    tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
                }
              for (unsigned int e = 0, k = dim; e < dim; ++e)
                for (unsigned int f = e + 1; f < dim; ++k, ++f)
                  for (unsigned int d = 0; d < dim; ++d)
                    {
                      tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
                      tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
                    }
              for (unsigned int d = 0; d < dim; ++d)
                {
                  tmp4[d] = tmp3[0] * inv_t_jac[d][0];
                  for (unsigned int e = 1; e < dim; ++e)
                    tmp4[d] += tmp3[e] * inv_t_jac[d][e];
                }
 
              VectorizedArrayType tmp2[dim];
              for (unsigned int d = 0; d < dim; ++d)
                {
                  tmp2[d] = t_jac[0][d] * this->values_quad[q_point];
                  for (unsigned e = 1; e < dim; ++e)
                    tmp2[d] +=
                      t_jac[e][d] * this->values_quad[e * nqp + q_point];
                }
 
              for (unsigned int d = 0; d < dim; ++d)
                for (unsigned int e = 0; e < dim; ++e)
                  {
                    grad_out[d][e] -= tmp4[e] * tmp2[d];
 
                    // finally multiply by det(J^{-1}) necessary in all
                    // contributions above
                    grad_out[d][e] *= inv_det;
                  }
            }
          return grad_out;
        }
    }
  Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_out;
 
  // Cartesian cell
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      for (unsigned int comp = 0; comp < n_components; ++comp)
        for (unsigned int d = 0; d < dim; ++d)
          grad_out[comp][d] =
            this->gradients_quad[(comp * nqp + q_point) * dim + d] *
            this->jacobian[0][d][d];
    }
  // cell with general/affine Jacobian
  else
    {
      const Tensor<2, dim, VectorizedArrayType> &jac =
        this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
                         q_point :
                         0];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        for (unsigned int d = 0; d < dim; ++d)
          {
            grad_out[comp][d] =
              jac[d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
            for (unsigned int e = 1; e < dim; ++e)
              grad_out[comp][d] +=
                jac[d][e] *
                this->gradients_quad[(comp * nqp + q_point) * dim + e];
          }
    }
  if constexpr (n_components == 1)
    return grad_out[0];
  else
    return grad_out;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE
  typename FEEvaluationBase<dim,
                            n_components_,
                            Number,
                            is_face,
                            VectorizedArrayType>::value_type
  FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
    get_normal_derivative(const unsigned int q_point) const
{
  AssertIndexRange(q_point, this->n_quadrature_points);
  if constexpr (running_in_debug_mode())
    {
      Assert(this->gradients_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
 
  Assert(this->normal_x_jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
 
  const std::size_t                            nqp = this->n_quadrature_points;
  Tensor<1, n_components, VectorizedArrayType> grad_out;
 
  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
    for (unsigned int comp = 0; comp < n_components; ++comp)
      grad_out[comp] =
        this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
        (this->normal_x_jacobian[0][dim - 1]);
  else
    {
      const std::size_t index =
        this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
                           this->normal_x_jacobian[index][0];
          for (unsigned int d = 1; d < dim; ++d)
            grad_out[comp] +=
              this->gradients_quad[(comp * nqp + q_point) * dim + d] *
              this->normal_x_jacobian[index][d];
        }
    }
  if constexpr (n_components == 1)
    return grad_out[0];
  else
    return grad_out;
}
 
 
 
namespace internal
{
  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
  // store the lower diagonal because of symmetry
  template <typename VectorizedArrayType>
  inline void
  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
                         const VectorizedArrayType *const         hessians,
                         const unsigned int,
                         VectorizedArrayType (&tmp)[1][1])
  {
    tmp[0][0] = jac[0][0] * hessians[0];
  }
 
  template <typename VectorizedArrayType>
  inline void
  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
                         const VectorizedArrayType *const         hessians,
                         const unsigned int                       nqp,
                         VectorizedArrayType (&tmp)[2][2])
  {
    for (unsigned int d = 0; d < 2; ++d)
      {
        tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
        tmp[1][d] =
          (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
      }
  }
 
  template <typename VectorizedArrayType>
  inline void
  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
                         const VectorizedArrayType *const         hessians,
                         const unsigned int                       nqp,
                         VectorizedArrayType (&tmp)[3][3])
  {
    for (unsigned int d = 0; d < 3; ++d)
      {
        tmp[0][d] =
          (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
           jac[d][2] * hessians[4 * nqp]);
        tmp[1][d] =
          (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
           jac[d][2] * hessians[5 * nqp]);
        tmp[2][d] =
          (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
           jac[d][2] * hessians[2 * nqp]);
      }
  }
} // namespace internal
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline typename FEEvaluationBase<dim,
                                 n_components_,
                                 Number,
                                 is_face,
                                 VectorizedArrayType>::hessian_type
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_hessian(const unsigned int q_point) const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->hessians_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
 
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_hessian"));
  const Tensor<2, dim, VectorizedArrayType> &jac =
    this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
                     0 :
                     q_point];
 
  Tensor<1, n_components, Tensor<2, dim, VectorizedArrayType>> hessian_out;
 
  const std::size_t      nqp  = this->n_quadrature_points;
  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
 
  // Cartesian cell
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          for (unsigned int d = 0; d < dim; ++d)
            hessian_out[comp][d][d] =
              this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
              (jac[d][d] * jac[d][d]);
          switch (dim)
            {
              case 1:
                break;
              case 2:
                hessian_out[comp][0][1] =
                  this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
                  (jac[0][0] * jac[1][1]);
                break;
              case 3:
                hessian_out[comp][0][1] =
                  this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
                  (jac[0][0] * jac[1][1]);
                hessian_out[comp][0][2] =
                  this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
                  (jac[0][0] * jac[2][2]);
                hessian_out[comp][1][2] =
                  this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
                  (jac[1][1] * jac[2][2]);
                break;
              default:
                DEAL_II_NOT_IMPLEMENTED();
            }
          for (unsigned int d = 0; d < dim; ++d)
            for (unsigned int e = d + 1; e < dim; ++e)
              hessian_out[comp][e][d] = hessian_out[comp][d][e];
        }
    }
  // cell with general Jacobian, but constant within the cell
  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          VectorizedArrayType tmp[dim][dim];
          internal::hessian_unit_times_jac(
            jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
 
          // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
          for (unsigned int d = 0; d < dim; ++d)
            for (unsigned int e = d; e < dim; ++e)
              {
                hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
                for (unsigned int f = 1; f < dim; ++f)
                  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
              }
 
          // no J' * grad(u) part here because the Jacobian is constant
          // throughout the cell and hence, its derivative is zero
 
          // take symmetric part
          for (unsigned int d = 0; d < dim; ++d)
            for (unsigned int e = d + 1; e < dim; ++e)
              hessian_out[comp][e][d] = hessian_out[comp][d][e];
        }
    }
  // cell with general Jacobian
  else
    {
      const auto &jac_grad = this->jacobian_gradients[q_point];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          VectorizedArrayType tmp[dim][dim];
          internal::hessian_unit_times_jac(
            jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
 
          // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
          for (unsigned int d = 0; d < dim; ++d)
            for (unsigned int e = d; e < dim; ++e)
              {
                hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
                for (unsigned int f = 1; f < dim; ++f)
                  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
              }
 
          // add diagonal part of J' * grad(u)
          for (unsigned int d = 0; d < dim; ++d)
            for (unsigned int e = 0; e < dim; ++e)
              hessian_out[comp][d][d] +=
                jac_grad[d][e] *
                this->gradients_quad[(comp * nqp + q_point) * dim + e];
 
          // add off-diagonal part of J' * grad(u)
          for (unsigned int d = 0, count = dim; d < dim; ++d)
            for (unsigned int e = d + 1; e < dim; ++e, ++count)
              for (unsigned int f = 0; f < dim; ++f)
                hessian_out[comp][d][e] +=
                  jac_grad[count][f] *
                  this->gradients_quad[(comp * nqp + q_point) * dim + f];
 
          // take symmetric part
          for (unsigned int d = 0; d < dim; ++d)
            for (unsigned int e = d + 1; e < dim; ++e)
              hessian_out[comp][e][d] = hessian_out[comp][d][e];
        }
    }
  if constexpr (n_components == 1)
    return hessian_out[0];
  else
    return hessian_out;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline typename FEEvaluationBase<dim,
                                 n_components_,
                                 Number,
                                 is_face,
                                 VectorizedArrayType>::gradient_type
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_hessian_diagonal(const unsigned int q_point) const
{
  Assert(!is_face, ExcNotImplemented());
  if constexpr (running_in_debug_mode())
    {
      Assert(this->hessians_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
 
  Assert(this->jacobian != nullptr, ExcNotImplemented());
  const Tensor<2, dim, VectorizedArrayType> &jac =
    this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
                     0 :
                     q_point];
 
  const std::size_t      nqp  = this->n_quadrature_points;
  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
  Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> hessian_out;
 
  // Cartesian cell
  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      for (unsigned int comp = 0; comp < n_components; ++comp)
        for (unsigned int d = 0; d < dim; ++d)
          hessian_out[comp][d] =
            this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
            (jac[d][d] * jac[d][d]);
    }
  // cell with general Jacobian, but constant within the cell
  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
    {
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          // compute laplacian before the gradient because it needs to access
          // unscaled gradient data
          VectorizedArrayType tmp[dim][dim];
          internal::hessian_unit_times_jac(
            jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
 
          // compute only the trace part of hessian, J * tmp = J *
          // hess_unit(u) * J^T
          for (unsigned int d = 0; d < dim; ++d)
            {
              hessian_out[comp][d] = jac[d][0] * tmp[0][d];
              for (unsigned int f = 1; f < dim; ++f)
                hessian_out[comp][d] += jac[d][f] * tmp[f][d];
            }
        }
    }
  // cell with general Jacobian
  else
    {
      const auto &jac_grad = this->jacobian_gradients[q_point];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          // compute laplacian before the gradient because it needs to access
          // unscaled gradient data
          VectorizedArrayType tmp[dim][dim];
          internal::hessian_unit_times_jac(
            jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
 
          // compute only the trace part of hessian, J * tmp = J *
          // hess_unit(u) * J^T
          for (unsigned int d = 0; d < dim; ++d)
            {
              hessian_out[comp][d] = jac[d][0] * tmp[0][d];
              for (unsigned int f = 1; f < dim; ++f)
                hessian_out[comp][d] += jac[d][f] * tmp[f][d];
            }
 
          for (unsigned int d = 0; d < dim; ++d)
            for (unsigned int e = 0; e < dim; ++e)
              hessian_out[comp][d] +=
                jac_grad[d][e] *
                this->gradients_quad[(comp * nqp + q_point) * dim + e];
        }
    }
 
  if constexpr (n_components == 1)
    return hessian_out[0];
  else
    return hessian_out;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline typename FEEvaluationBase<dim,
                                 n_components_,
                                 Number,
                                 is_face,
                                 VectorizedArrayType>::value_type
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_laplacian(const unsigned int q_point) const
{
  Assert(is_face == false, ExcNotImplemented());
  if constexpr (running_in_debug_mode())
    {
      Assert(this->hessians_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
 
  const gradient_type hess_diag = get_hessian_diagonal(q_point);
  if constexpr (n_components == 1)
    {
      VectorizedArrayType sum = hess_diag[0];
      for (unsigned int d = 1; d < dim; ++d)
        sum += hess_diag[d];
      return sum;
    }
  else
    {
      Tensor<1, n_components_, VectorizedArrayType> laplacian_out;
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          laplacian_out[comp] = hess_diag[comp][0];
          for (unsigned int d = 1; d < dim; ++d)
            laplacian_out[comp] += hess_diag[comp][d];
        }
      return laplacian_out;
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline typename FEEvaluationBase<dim,
                                 n_components_,
                                 Number,
                                 is_face,
                                 VectorizedArrayType>::value_type
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_normal_hessian(const unsigned int q_point) const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->hessians_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
 
  Assert(this->normal_x_jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_hessians"));
 
  Tensor<1, n_components, VectorizedArrayType> hessian_out;
 
  const std::size_t      nqp  = this->n_quadrature_points;
  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
 
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      const auto nxj = this->normal_x_jacobian[0];
 
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          for (unsigned int d = 0; d < dim; ++d)
            hessian_out[comp] +=
              this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
              (nxj[d]) * (nxj[d]);
 
          switch (dim)
            {
              case 1:
                break;
              case 2:
                hessian_out[comp] +=
                  this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
                  (nxj[0] * nxj[1]);
                break;
              case 3:
                hessian_out[comp] +=
                  2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
                  (nxj[0] * nxj[1]);
                hessian_out[comp] +=
                  2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
                  (nxj[0] * nxj[2]);
                hessian_out[comp] +=
                  2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
                  (nxj[1] * nxj[2]);
                break;
              default:
                DEAL_II_NOT_IMPLEMENTED();
            }
        }
    }
  // cell with general Jacobian
  else
    {
      const auto normal  = this->normal_vector(q_point);
      const auto hessian = get_hessian(q_point);
 
      if constexpr (n_components == 1)
        hessian_out[0] = hessian * normal * normal;
      else
        for (unsigned int comp = 0; comp < n_components; ++comp)
          hessian_out[comp] = hessian[comp] * normal * normal;
    }
  if constexpr (n_components == 1)
    return hessian_out[0];
  else
    return hessian_out;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_dof_value(const value_type val_in, const unsigned int dof)
{
  if constexpr (running_in_debug_mode())
    {
      this->dof_values_initialized = true;
    }
  const std::size_t dofs = this->data->dofs_per_component_on_cell;
  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
  for (unsigned int comp = 0; comp < n_components; ++comp)
    if constexpr (n_components == 1)
      this->values_dofs[comp * dofs + dof] = val_in;
    else
      this->values_dofs[comp * dofs + dof] = val_in[comp];
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_value(const value_type val_in, const unsigned int q_point)
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->is_reinitialized, ExcNotInitialized());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_values"));
  if constexpr (running_in_debug_mode())
    {
      this->values_quad_submitted = true;
    }
 
  const std::size_t    nqp    = this->n_quadrature_points;
  VectorizedArrayType *values = this->values_quad + q_point;
 
  const VectorizedArrayType JxW =
    this->cell_type <= internal::MatrixFreeFunctions::affine ?
      this->J_value[0] * this->quadrature_weights[q_point] :
      this->J_value[q_point];
  if constexpr (n_components == 1)
    values[0] = val_in * JxW;
  else
    {
      if (n_components == dim &&
          this->data->element_type ==
            internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
        {
          // Piola transform is required
          AssertIndexRange(q_point, this->n_quadrature_points);
          Assert(this->J_value != nullptr,
                 internal::ExcMatrixFreeAccessToUninitializedMappingField(
                   "update_value"));
          if constexpr (running_in_debug_mode())
            {
              Assert(this->is_reinitialized, ExcNotInitialized());
              this->values_quad_submitted = true;
            }
 
          VectorizedArrayType *values = this->values_quad + q_point;
          const std::size_t    nqp    = this->n_quadrature_points;
 
          if (!is_face &&
              this->cell_type == internal::MatrixFreeFunctions::cartesian)
            {
              const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
              const VectorizedArrayType                 weight =
                this->quadrature_weights[q_point];
 
              for (unsigned int comp = 0; comp < n_components; ++comp)
                values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
            }
          else
            {
              // Affine or general cell
              const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
                (this->cell_type > internal::MatrixFreeFunctions::affine) ?
                  this->jacobian[q_point] :
                  this->jacobian[0];
 
              // Derivatives are reordered for faces. Need to take this into
              // account and 1/inv_det != J_value for faces
              const VectorizedArrayType fac =
                (!is_face) ?
                  this->quadrature_weights[q_point] :
                  (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
                      this->J_value[q_point] :
                      this->J_value[0] * this->quadrature_weights[q_point]) *
                   ((dim == 2 && this->get_face_no() < 2) ?
                      -determinant(inv_t_jac) :
                      determinant(inv_t_jac)));
              const Tensor<2, dim, VectorizedArrayType> jac =
                (this->cell_type > internal::MatrixFreeFunctions::affine) ?
                  transpose(invert(inv_t_jac)) :
                  this->jacobian[1];
 
              // J^T * u * factor
              for (unsigned int comp = 0; comp < n_components; ++comp)
                {
                  values[comp * nqp] = val_in[0] * jac[0][comp];
                  for (unsigned int e = 1; e < dim; ++e)
                    values[comp * nqp] += val_in[e] * jac[e][comp];
                  values[comp * nqp] *= fac;
                }
            }
        }
      else if (n_components == dim &&
               this->data->element_type ==
                 internal::MatrixFreeFunctions::ElementType::tensor_nedelec)
        {
          // Piola transform is required
          AssertIndexRange(q_point, this->n_quadrature_points);
          Assert(this->J_value != nullptr,
                 internal::ExcMatrixFreeAccessToUninitializedMappingField(
                   "update_value"));
          if constexpr (running_in_debug_mode())
            {
              Assert(this->is_reinitialized, ExcNotInitialized());
              this->values_quad_submitted = true;
            }
 
          VectorizedArrayType *values = this->values_quad + q_point;
          const std::size_t    nqp    = this->n_quadrature_points;
 
          if (!is_face &&
              this->cell_type == internal::MatrixFreeFunctions::cartesian)
            {
              const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
                this->jacobian[0];
              const VectorizedArrayType det_jac =
                (dim == 2) ? this->jacobian[1][0][0] * this->jacobian[1][1][1] :
                             this->jacobian[1][0][0] * this->jacobian[1][1][1] *
                               this->jacobian[1][2][2];
 
              const VectorizedArrayType weight =
                this->quadrature_weights[q_point];
 
              // J^{-1} * u * det(J) * weight
              for (unsigned int comp = 0; comp < n_components; ++comp)
                values[comp * nqp] =
                  val_in[comp] * weight * inv_t_jac[comp][comp] * det_jac;
            }
          else
            {
              // Affine or general cell
              const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
                (this->cell_type > internal::MatrixFreeFunctions::affine) ?
                  this->jacobian[q_point] :
                  this->jacobian[0];
 
              const VectorizedArrayType fac =
                ((this->cell_type > internal::MatrixFreeFunctions::affine) ?
                   this->J_value[q_point] :
                   this->J_value[0] * this->quadrature_weights[q_point]);
 
              // J^{-1} * u * factor
              for (unsigned int comp = 0; comp < n_components; ++comp)
                {
                  values[comp * nqp] = val_in[0] * inv_t_jac[0][comp];
                  for (unsigned int e = 1; e < dim; ++e)
                    values[comp * nqp] += val_in[e] * inv_t_jac[e][comp];
                  values[comp * nqp] *= fac;
                }
            }
        }
      else
        for (unsigned int comp = 0; comp < n_components; ++comp)
          values[comp * nqp] = val_in[comp] * JxW;
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <int, typename>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_value(const Tensor<1, 1, VectorizedArrayType> val_in,
               const unsigned int                      q_point)
{
  static_assert(n_components == 1,
                "Do not try to modify the default template parameters used for"
                " selectively enabling this function via std::enable_if!");
  submit_value(val_in[0], q_point);
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_gradient(const gradient_type grad_in, const unsigned int q_point)
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->is_reinitialized, ExcNotInitialized());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  if constexpr (running_in_debug_mode())
    {
      this->gradients_quad_submitted = true;
    }
 
  if constexpr (dim > 1 && n_components == dim)
    {
      if (this->data->element_type ==
          internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
        {
          // Piola transform is required
 
          if constexpr (running_in_debug_mode())
            {
              Assert(this->is_reinitialized, ExcNotInitialized());
            }
          AssertIndexRange(q_point, this->n_quadrature_points);
          Assert(this->J_value != nullptr,
                 internal::ExcMatrixFreeAccessToUninitializedMappingField(
                   "update_gradients"));
          Assert(this->jacobian != nullptr,
                 internal::ExcMatrixFreeAccessToUninitializedMappingField(
                   "update_gradients"));
          if constexpr (running_in_debug_mode())
            {
              this->gradients_quad_submitted = true;
            }
 
          VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
          VectorizedArrayType *values =
            this->values_from_gradients_quad + q_point;
          const std::size_t nqp   = this->n_quadrature_points;
          const std::size_t nqp_d = nqp * dim;
 
          if (!is_face &&
              this->cell_type == internal::MatrixFreeFunctions::cartesian)
            {
              // Cartesian cell
              const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
                this->jacobian[0];
              const Tensor<2, dim, VectorizedArrayType> &jac =
                this->jacobian[1];
              const VectorizedArrayType weight =
                this->quadrature_weights[q_point];
              for (unsigned int d = 0; d < dim; ++d)
                for (unsigned int comp = 0; comp < n_components; ++comp)
                  gradients[comp * nqp_d + d] = grad_in[comp][d] *
                                                inv_t_jac[d][d] *
                                                (jac[comp][comp] * weight);
            }
          else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
            {
              // Affine cell
              const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
                this->jacobian[0];
              const Tensor<2, dim, VectorizedArrayType> &jac =
                this->jacobian[1];
 
              // Derivatives are reordered for faces. Need to take this into
              // account and 1/inv_det != J_value for faces
              const VectorizedArrayType fac =
                (!is_face) ?
                  this->quadrature_weights[q_point] :
                  this->J_value[0] * this->quadrature_weights[q_point] *
                    ((dim == 2 && this->get_face_no() < 2) ?
                       -determinant(inv_t_jac) :
                       determinant(inv_t_jac));
 
              // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
              VectorizedArrayType tmp[dim][dim];
              for (unsigned int d = 0; d < dim; ++d)
                for (unsigned int e = 0; e < dim; ++e)
                  {
                    tmp[d][e] = inv_t_jac[0][d] * grad_in[e][0];
                    for (unsigned int f = 1; f < dim; ++f)
                      tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
                  }
              for (unsigned int comp = 0; comp < n_components; ++comp)
                for (unsigned int d = 0; d < dim; ++d)
                  {
                    VectorizedArrayType res = jac[0][comp] * tmp[d][0];
                    for (unsigned int f = 1; f < dim; ++f)
                      res += jac[f][comp] * tmp[d][f];
 
                    gradients[comp * nqp_d + d] = res * fac;
                  }
            }
          else
            {
              // General cell
 
              const auto jac_grad =
                this->jacobian_gradients_non_inverse[q_point];
              const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
                this->jacobian[q_point];
 
              // Derivatives are reordered for faces. Need to take this into
              // account and 1/inv_det != J_value for faces
              const VectorizedArrayType fac =
                (!is_face) ? this->quadrature_weights[q_point] :
                             this->J_value[q_point] *
                               ((dim == 2 && this->get_face_no() < 2) ?
                                  -determinant(inv_t_jac) :
                                  determinant(inv_t_jac));
              const Tensor<2, dim, VectorizedArrayType> t_jac =
                invert(inv_t_jac);
 
              // Start evaluation for values part below to enable the compiler
              // to possibly re-use the same computation in get_gradient()
              // without interfering with stores to 'gradients'
              VectorizedArrayType tmp3[dim], tmp4[dim];
              for (unsigned int d = 0; d < dim; ++d)
                {
                  tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
                  for (unsigned int e = 1; e < dim; ++e)
                    tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
                }
              for (unsigned int e = 0, k = dim; e < dim; ++e)
                for (unsigned int f = e + 1; f < dim; ++k, ++f)
                  for (unsigned int d = 0; d < dim; ++d)
                    {
                      tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
                      tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
                    }
              for (unsigned int d = 0; d < dim; ++d)
                {
                  tmp4[d] = tmp3[0] * inv_t_jac[d][0];
                  for (unsigned int e = 1; e < dim; ++e)
                    tmp4[d] += tmp3[e] * inv_t_jac[d][e];
                }
 
              const Tensor<2, dim, VectorizedArrayType> grad_in_scaled =
                fac * grad_in;
 
              VectorizedArrayType tmp[dim][dim];
 
              // J * (J^{-1} * (grad_in * factor))
              for (unsigned int d = 0; d < dim; ++d)
                for (unsigned int e = 0; e < dim; ++e)
                  {
                    tmp[d][e] = inv_t_jac[0][d] * grad_in_scaled[e][0];
                    for (unsigned int f = 1; f < dim; ++f)
                      tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
                  }
 
              for (unsigned int d = 0; d < dim; ++d)
                for (unsigned int e = 0; e < dim; ++e)
                  {
                    VectorizedArrayType res = t_jac[d][0] * tmp[e][0];
                    for (unsigned int f = 1; f < dim; ++f)
                      res += t_jac[d][f] * tmp[e][f];
 
                    gradients[d * nqp_d + e] = res;
                  }
 
              // jac_grad * (J^{-1} * (grad_in * factor)), re-use part in braces
              // as 'tmp' from above
              VectorizedArrayType value[dim];
              for (unsigned int d = 0; d < dim; ++d)
                {
                  value[d] = tmp[d][0] * jac_grad[d][0];
                  for (unsigned int e = 1; e < dim; ++e)
                    value[d] += tmp[d][e] * jac_grad[d][e];
                }
              for (unsigned int e = 0, k = dim; e < dim; ++e)
                for (unsigned int f = e + 1; f < dim; ++k, ++f)
                  for (unsigned int d = 0; d < dim; ++d)
                    {
                      value[e] += tmp[f][d] * jac_grad[k][d];
                      value[f] += tmp[e][d] * jac_grad[k][d];
                    }
 
              //   -(grad_in * factor) * J * (J^{-T} * jac_grad * J^{-1})
              // = -(grad_in * factor) * J * ( \------- tmp4 ---------/ )
              for (unsigned int d = 0; d < dim; ++d)
                {
                  VectorizedArrayType tmp2 = grad_in_scaled[d][0] * tmp4[0];
                  for (unsigned int e = 1; e < dim; ++e)
                    tmp2 += grad_in_scaled[d][e] * tmp4[e];
                  for (unsigned int e = 0; e < dim; ++e)
                    value[e] -= t_jac[e][d] * tmp2;
                }
 
              for (unsigned int d = 0; d < dim; ++d)
                values[d * nqp] = value[d];
            }
          return;
        }
    }
 
  const std::size_t    nqp_d     = this->n_quadrature_points * dim;
  VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
 
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
 
      // Load all entries before starting to write back to make sure the
      // compiler sees opportunity of loads in a possibly nearby
      // get_gradient() function (i.e., the compiler should not think that
      // 'jacobian' could alias with 'gradients').
      std::array<VectorizedArrayType, dim> jac;
      for (unsigned int d = 0; d < dim; ++d)
        jac[d] = this->jacobian[0][d][d];
 
      for (unsigned int d = 0; d < dim; ++d)
        {
          const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
          if constexpr (n_components == 1)
            gradients[d] = grad_in[d] * factor;
          else
            for (unsigned int comp = 0; comp < n_components; ++comp)
              gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
        }
    }
  else
    {
      const Tensor<2, dim, VectorizedArrayType> jac =
        this->cell_type > internal::MatrixFreeFunctions::affine ?
          this->jacobian[q_point] :
          this->jacobian[0];
      const VectorizedArrayType JxW =
        this->cell_type > internal::MatrixFreeFunctions::affine ?
          this->J_value[q_point] :
          this->J_value[0] * this->quadrature_weights[q_point];
      if constexpr (n_components == 1)
        for (unsigned int d = 0; d < dim; ++d)
          {
            VectorizedArrayType new_val = jac[0][d] * grad_in[0];
            for (unsigned int e = 1; e < dim; ++e)
              new_val += (jac[e][d] * grad_in[e]);
            gradients[d] = new_val * JxW;
          }
      else
        for (unsigned int comp = 0; comp < n_components; ++comp)
          for (unsigned int d = 0; d < dim; ++d)
            {
              VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
              for (unsigned int e = 1; e < dim; ++e)
                new_val += (jac[e][d] * grad_in[comp][e]);
              gradients[comp * nqp_d + d] = new_val * JxW;
            }
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <int, typename>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_gradient(const Tensor<2, 1, VectorizedArrayType> grad_in,
                  const unsigned int                      q_point)
{
  static_assert(n_components == 1 && dim == 1,
                "Do not try to modify the default template parameters used for"
                " selectively enabling this function via std::enable_if!");
  submit_gradient(grad_in[0], q_point);
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
{
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->normal_x_jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  if constexpr (running_in_debug_mode())
    {
      this->gradients_quad_submitted = true;
    }
 
  const std::size_t    nqp_d     = this->n_quadrature_points * dim;
  VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
 
  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArrayType JxW_jac = this->J_value[0] *
                                          this->quadrature_weights[q_point] *
                                          this->normal_x_jacobian[0][dim - 1];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          for (unsigned int d = 0; d < dim - 1; ++d)
            gradients[comp * nqp_d + d] = VectorizedArrayType();
          if constexpr (n_components == 1)
            gradients[dim - 1] = grad_in * JxW_jac;
          else
            gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
        }
    }
  else
    {
      const unsigned int index =
        this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
      const Tensor<1, dim, VectorizedArrayType> jac =
        this->normal_x_jacobian[index];
      const VectorizedArrayType JxW =
        (this->cell_type <= internal::MatrixFreeFunctions::affine) ?
          this->J_value[index] * this->quadrature_weights[q_point] :
          this->J_value[index];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        for (unsigned int d = 0; d < dim; ++d)
          if constexpr (n_components == 1)
            gradients[d] = (grad_in * JxW) * jac[d];
          else
            gradients[comp * nqp_d + d] = (grad_in[comp] * JxW) * jac[d];
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->is_reinitialized, ExcNotInitialized());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_hessians"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_hessians"));
  if constexpr (running_in_debug_mode())
    {
      this->hessians_quad_submitted = true;
    }
 
  // compute hessian_unit = J^T * hessian_in(u) * J
  const std::size_t      nqp  = this->n_quadrature_points;
  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
 
      // diagonal part
      for (unsigned int d = 0; d < dim; ++d)
        {
          const auto                jac_d  = this->jacobian[0][d][d];
          const VectorizedArrayType factor = jac_d * jac_d * JxW;
          for (unsigned int comp = 0; comp < n_components; ++comp)
            if constexpr (n_components == 1)
              this->hessians_quad[d * nqp + q_point] =
                hessian_in[d][d] * factor;
            else
              this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
                hessian_in[comp][d][d] * factor;
        }
 
      // off diagonal part
      for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
        for (unsigned int e = 0; e < d; ++e, ++off_dia)
          {
            const auto                jac_d  = this->jacobian[0][d][d];
            const auto                jac_e  = this->jacobian[0][e][e];
            const VectorizedArrayType factor = jac_d * jac_e * JxW;
            for (unsigned int comp = 0; comp < n_components; ++comp)
              if constexpr (n_components == 1)
                this->hessians_quad[off_dia * nqp + q_point] =
                  (hessian_in[d][e] + hessian_in[e][d]) * factor;
              else
                this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
                  (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
          }
    }
  // cell with general Jacobian, but constant within the cell
  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
      const VectorizedArrayType                 JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          Tensor<2, dim, VectorizedArrayType> hessian_c;
          if constexpr (n_components == 1)
            hessian_c = hessian_in;
          else
            hessian_c = hessian_in[comp];
 
          // 1. tmp = hessian(u) * J
          VectorizedArrayType tmp[dim][dim];
          for (unsigned int i = 0; i < dim; ++i)
            for (unsigned int j = 0; j < dim; ++j)
              {
                tmp[i][j] = hessian_c[i][0] * jac[0][j];
                for (unsigned int k = 1; k < dim; ++k)
                  tmp[i][j] += hessian_c[i][k] * jac[k][j];
              }
 
          // 2. hessian_unit = J^T * tmp
          VectorizedArrayType tmp2[dim][dim];
          for (unsigned int i = 0; i < dim; ++i)
            for (unsigned int j = 0; j < dim; ++j)
              {
                tmp2[i][j] = jac[0][i] * tmp[0][j];
                for (unsigned int k = 1; k < dim; ++k)
                  tmp2[i][j] += jac[k][i] * tmp[k][j];
              }
 
          // diagonal part
          for (unsigned int d = 0; d < dim; ++d)
            this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
              tmp2[d][d] * JxW;
 
          // off diagonal part
          for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
            for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
              this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
                (tmp2[d][e] + tmp2[e][d]) * JxW;
        }
    }
  else
    {
      const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
      const VectorizedArrayType                 JxW = this->J_value[q_point];
      const auto &jac_grad = this->jacobian_gradients[q_point];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        {
          Tensor<2, dim, VectorizedArrayType> hessian_c;
          if constexpr (n_components == 1)
            hessian_c = hessian_in;
          else
            hessian_c = hessian_in[comp];
 
          // 1. tmp = hessian(u) * J
          VectorizedArrayType tmp[dim][dim];
          for (unsigned int i = 0; i < dim; ++i)
            for (unsigned int j = 0; j < dim; ++j)
              {
                tmp[i][j] = hessian_c[i][0] * jac[0][j];
                for (unsigned int k = 1; k < dim; ++k)
                  tmp[i][j] += hessian_c[i][k] * jac[k][j];
              }
 
          // 2. hessian_unit = J^T * tmp
          VectorizedArrayType tmp2[dim][dim];
          for (unsigned int i = 0; i < dim; ++i)
            for (unsigned int j = 0; j < dim; ++j)
              {
                tmp2[i][j] = jac[0][i] * tmp[0][j];
                for (unsigned int k = 1; k < dim; ++k)
                  tmp2[i][j] += jac[k][i] * tmp[k][j];
              }
 
          // diagonal part
          for (unsigned int d = 0; d < dim; ++d)
            this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
              tmp2[d][d] * JxW;
 
          // off diagonal part
          for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
            for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
              this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
                (tmp2[d][e] + tmp2[e][d]) * JxW;
 
          // 3. gradient_unit = J' * hessian
          for (unsigned int d = 0; d < dim; ++d)
            {
              VectorizedArrayType sum = 0;
              for (unsigned int e = 0; e < dim; ++e)
                sum += hessian_c[e][e] * jac_grad[e][d];
              for (unsigned int e = 0, count = dim; e < dim; ++e)
                for (unsigned int f = e + 1; f < dim; ++f, ++count)
                  sum +=
                    (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][d];
              this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
                                                 d] = sum * JxW;
            }
        }
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_normal_hessian(const value_type   normal_hessian_in,
                        const unsigned int q_point)
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->is_reinitialized, ExcNotInitialized());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_hessians"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_hessians"));
  if constexpr (running_in_debug_mode())
    {
      this->hessians_quad_submitted = true;
    }
 
  // compute hessian_unit = J^T * hessian_in(u) * J
  const std::size_t      nqp  = this->n_quadrature_points;
  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
 
      const auto nxj = this->normal_x_jacobian[0];
 
      // diagonal part
      for (unsigned int d = 0; d < dim; ++d)
        {
          const auto                nxj_d  = nxj[d];
          const VectorizedArrayType factor = nxj_d * nxj_d * JxW;
          for (unsigned int comp = 0; comp < n_components; ++comp)
            if constexpr (n_components == 1)
              this->hessians_quad[d * nqp + q_point] =
                normal_hessian_in * factor;
            else
              this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
                normal_hessian_in[comp] * factor;
        }
 
      // off diagonal part
      for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
        for (unsigned int e = 0; e < d; ++e, ++off_dia)
          {
            const auto                jac_d  = nxj[d];
            const auto                jac_e  = nxj[e];
            const VectorizedArrayType factor = jac_d * jac_e * JxW;
            for (unsigned int comp = 0; comp < n_components; ++comp)
              if constexpr (n_components == 1)
                this->hessians_quad[off_dia * nqp + q_point] =
                  2. * normal_hessian_in * factor;
              else
                this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
                  2. * normal_hessian_in[comp] * factor;
          }
    }
  else
    {
      const auto normal           = this->normal_vector(q_point);
      const auto normal_projector = outer_product(normal, normal);
      if constexpr (n_components == 1)
        submit_hessian(normal_hessian_in * normal_projector, q_point);
      else
        {
          hessian_type tmp;
          for (unsigned int comp = 0; comp < n_components; ++comp)
            tmp[comp] = normal_hessian_in[comp] * normal_projector;
          submit_hessian(tmp, q_point);
        }
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline typename FEEvaluationBase<dim,
                                 n_components_,
                                 Number,
                                 is_face,
                                 VectorizedArrayType>::value_type
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  integrate_value() const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->is_reinitialized, ExcNotInitialized());
      Assert(this->values_quad_submitted == true,
             internal::ExcAccessToUninitializedField());
    }
 
  Tensor<1, n_components, VectorizedArrayType> return_value;
  const std::size_t                            nqp = this->n_quadrature_points;
  for (unsigned int q = 0; q < nqp; ++q)
    for (unsigned int comp = 0; comp < n_components; ++comp)
      return_value[comp] += this->values_quad[comp * nqp + q];
  if constexpr (n_components == 1)
    return return_value[0];
  else
    return return_value;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <int, typename>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_divergence(const unsigned int q_point) const
{
  static_assert(n_components == dim,
                "Do not try to modify the default template parameters used for"
                " selectively enabling this function via std::enable_if!");
 
  if constexpr (running_in_debug_mode())
    {
      Assert(this->gradients_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
 
  VectorizedArrayType divergence;
  const std::size_t   nqp = this->n_quadrature_points;
 
  if (dim > 1 &&
      this->data->element_type ==
        internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
    {
      VectorizedArrayType inv_det =
        (!is_face &&
         this->cell_type == internal::MatrixFreeFunctions::cartesian) ?
          this->jacobian[0][0][0] *
            ((dim == 2) ? this->jacobian[0][1][1] :
                          this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
          determinant(this->jacobian[this->cell_type >
                                         internal::MatrixFreeFunctions::affine ?
                                       q_point :
                                       0]);
 
      // on faces in 2d, the determinant has the wrong sign due to ordering of
      // derivatives
      if (is_face && dim == 2 && this->get_face_no() < 2)
        inv_det = -inv_det;
 
      // div * det(J^-1)
      divergence = this->gradients_quad[q_point * dim];
      for (unsigned int d = 1; d < dim; ++d)
        divergence += this->gradients_quad[(d * nqp + q_point) * dim + d];
      divergence *= inv_det;
    }
  else
    {
      if (!is_face &&
          this->cell_type == internal::MatrixFreeFunctions::cartesian)
        {
          // Cartesian cell
          divergence =
            this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
          for (unsigned int d = 1; d < dim; ++d)
            divergence += this->gradients_quad[(d * nqp + q_point) * dim + d] *
                          this->jacobian[0][d][d];
        }
      else
        {
          // cell with general/constant Jacobian
          const Tensor<2, dim, VectorizedArrayType> &jac =
            this->cell_type == internal::MatrixFreeFunctions::general ?
              this->jacobian[q_point] :
              this->jacobian[0];
          divergence = jac[0][0] * this->gradients_quad[q_point * dim];
          for (unsigned int e = 1; e < dim; ++e)
            divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
          for (unsigned int d = 1; d < dim; ++d)
            for (unsigned int e = 0; e < dim; ++e)
              divergence +=
                jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim + e];
        }
    }
  return divergence;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <int, typename>
inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, VectorizedArrayType>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_symmetric_gradient(const unsigned int q_point) const
{
  static_assert(n_components == dim,
                "Do not try to modify the default template parameters used for"
                " selectively enabling this function via std::enable_if!");
 
  // copy from generic function into dim-specialization function
  const auto          grad = get_gradient(q_point);
  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
  VectorizedArrayType half = Number(0.5);
  for (unsigned int d = 0; d < dim; ++d)
    symmetrized[d] = grad[d][d];
  switch (dim)
    {
      case 1:
        break;
      case 2:
        symmetrized[2] = grad[0][1] + grad[1][0];
        symmetrized[2] *= half;
        break;
      case 3:
        symmetrized[3] = grad[0][1] + grad[1][0];
        symmetrized[3] *= half;
        symmetrized[4] = grad[0][2] + grad[2][0];
        symmetrized[4] *= half;
        symmetrized[5] = grad[1][2] + grad[2][1];
        symmetrized[5] *= half;
        break;
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
  return SymmetricTensor<2, dim, VectorizedArrayType>(symmetrized);
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <int, typename>
inline DEAL_II_ALWAYS_INLINE
  typename FEEvaluationBase<dim,
                            n_components_,
                            Number,
                            is_face,
                            VectorizedArrayType>::curl_type
  FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
    get_curl(const unsigned int q_point) const
{
  static_assert(dim > 1 && n_components == dim,
                "Do not try to modify the default template parameters used for"
                " selectively enabling this function via std::enable_if!");
 
  if (dim > 1 && this->data->element_type ==
                   internal::MatrixFreeFunctions::ElementType::tensor_nedelec)
    {
      // Piola transformation required: Since the transformation for Nedelec
      // is constructed to give a curl based on a scaled version of the
      // reference curl, this function turns out to be particularly simple.
      if constexpr (running_in_debug_mode())
        {
          Assert(this->gradients_quad_initialized == true,
                 internal::ExcAccessToUninitializedField());
        }
      AssertIndexRange(q_point, this->n_quadrature_points);
      Assert(this->jacobian != nullptr,
             internal::ExcMatrixFreeAccessToUninitializedMappingField(
               "update_gradients"));
      curl_type                  curl;
      const std::size_t          nqp   = this->n_quadrature_points;
      const std::size_t          nqp_d = nqp * dim;
      const VectorizedArrayType *gradients =
        this->gradients_quad + q_point * dim;
      if constexpr (dim == 2)
        curl = gradients[1 * nqp_d + 0] - gradients[0 * nqp_d + 1];
      else if constexpr (dim == 3)
        {
          curl[0] = gradients[2 * nqp_d + 1] - gradients[1 * nqp_d + 2];
          curl[1] = gradients[0 * nqp_d + 2] - gradients[2 * nqp_d + 0];
          curl[2] = gradients[1 * nqp_d + 0] - gradients[0 * nqp_d + 1];
        }
      else
        DEAL_II_NOT_IMPLEMENTED();
 
      if (!is_face &&
          this->cell_type == internal::MatrixFreeFunctions::cartesian)
        {
          VectorizedArrayType inv_det =
            this->jacobian[0][0][0] *
            ((dim == 2) ? this->jacobian[0][1][1] :
                          this->jacobian[0][1][1] * this->jacobian[0][2][2]);
 
          const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
 
          // J * curl * det(J^-1)
          if constexpr (dim == 2)
            curl *= inv_det;
          else
            for (unsigned int d = 0; d < (dim == 2 ? 1 : dim); ++d)
              curl[d] *= jac[d][d] * inv_det;
        }
      else
        {
          // Cell with general/constant Jacobian.
          const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
            (this->cell_type > internal::MatrixFreeFunctions::affine) ?
              this->jacobian[q_point] :
              this->jacobian[0];
          const Tensor<2, dim, VectorizedArrayType> jac =
            (this->cell_type > internal::MatrixFreeFunctions::affine) ?
              transpose(invert(inv_t_jac)) :
              this->jacobian[1];
 
          VectorizedArrayType inv_det = determinant(inv_t_jac);
 
          // on faces in 2d, the determinant has the wrong sign due to ordering
          // of derivatives
          if (is_face && dim == 2 && this->get_face_no() < 2)
            inv_det = -inv_det;
 
          // J * curl * det(J^-1)
          // In case of 2D its just a scaling with det(J^-1)
          if constexpr (dim == 2)
            curl *= inv_det;
          else
            {
              Tensor<1, dim, VectorizedArrayType> curl_temp;
              for (unsigned int comp = 0; comp < dim; ++comp)
                {
                  curl_temp[comp] = curl[0] * jac[comp][0];
                  for (unsigned int e = 1; e < dim; ++e)
                    curl_temp[comp] += jac[comp][e] * curl[e];
                }
              for (unsigned int comp = 0; comp < dim; ++comp)
                curl[comp] = curl_temp[comp] * inv_det;
            }
        }
      return curl;
    }
 
  // copy from generic function into dim-specialization function
  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
  curl_type                                 curl;
  if constexpr (dim == 2)
    curl = grad[1][0] - grad[0][1];
  else if constexpr (dim == 3)
    {
      curl[0] = grad[2][1] - grad[1][2];
      curl[1] = grad[0][2] - grad[2][0];
      curl[2] = grad[1][0] - grad[0][1];
    }
  else
    DEAL_II_NOT_IMPLEMENTED();
 
  return curl;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <int, typename>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_divergence(const VectorizedArrayType div_in,
                    const unsigned int        q_point)
{
  static_assert(n_components == dim,
                "Do not try to modify the default template parameters used for"
                " selectively enabling this function via std::enable_if!");
 
  if constexpr (running_in_debug_mode())
    {
      Assert(this->is_reinitialized, ExcNotInitialized());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  if constexpr (running_in_debug_mode())
    {
      this->gradients_quad_submitted = true;
    }
 
  const std::size_t    nqp_d     = this->n_quadrature_points * dim;
  VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
 
  if (this->data->element_type ==
      internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
    {
      // General cell
 
      // Derivatives are reordered for faces. Need to take this into account
      // and 1/inv_det != J_value for faces
      const VectorizedArrayType fac =
        (!is_face) ?
          this->quadrature_weights[q_point] * div_in :
          (this->cell_type > internal::MatrixFreeFunctions::affine ?
             this->J_value[q_point] :
             this->J_value[0] * this->quadrature_weights[q_point]) *
            div_in *
            determinant(
              this->jacobian[this->cell_type >
                                 internal::MatrixFreeFunctions::affine ?
                               q_point :
                               0]) *
            Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
 
      for (unsigned int d = 0; d < dim; ++d)
        {
          for (unsigned int e = 0; e < dim; ++e)
            gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
        }
      this->divergence_is_requested = true;
    }
  else
    {
      if (!is_face &&
          this->cell_type == internal::MatrixFreeFunctions::cartesian)
        {
          const VectorizedArrayType fac =
            this->J_value[0] * this->quadrature_weights[q_point] * div_in;
          for (unsigned int d = 0; d < dim; ++d)
            {
              const VectorizedArrayType jac_dd = this->jacobian[0][d][d];
              for (unsigned int e = 0; e < dim; ++e)
                gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
            }
        }
      else
        {
          const Tensor<2, dim, VectorizedArrayType> jac =
            this->cell_type == internal::MatrixFreeFunctions::general ?
              this->jacobian[q_point] :
              this->jacobian[0];
          const VectorizedArrayType fac =
            (this->cell_type == internal::MatrixFreeFunctions::general ?
               this->J_value[q_point] :
               this->J_value[0] * this->quadrature_weights[q_point]) *
            div_in;
          for (unsigned int d = 0; d < dim; ++d)
            {
              for (unsigned int e = 0; e < dim; ++e)
                gradients[d * nqp_d + e] = jac[d][e] * fac;
            }
        }
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <int, typename>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_symmetric_gradient(
    const SymmetricTensor<2, dim, VectorizedArrayType> sym_grad,
    const unsigned int                                 q_point)
{
  static_assert(n_components == dim,
                "Do not try to modify the default template parameters used for"
                " selectively enabling this function via std::enable_if!");
 
  AssertThrow(
    this->data->element_type !=
      internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas,
    ExcNotImplemented());
 
  // could have used base class operator, but that involves some overhead
  // which is inefficient. it is nice to have the symmetric tensor because
  // that saves some operations
  if constexpr (running_in_debug_mode())
    {
      Assert(this->is_reinitialized, ExcNotInitialized());
    }
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  if constexpr (running_in_debug_mode())
    {
      this->gradients_quad_submitted = true;
    }
 
  const std::size_t    nqp_d     = this->n_quadrature_points * dim;
  VectorizedArrayType *gradients = this->gradients_quad + dim * q_point;
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
      const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
      for (unsigned int d = 0; d < dim; ++d)
        gradients[d * nqp_d + d] =
          (sym_grad.access_raw_entry(d) * JxW * jac[d][d]);
      for (unsigned int e = 0, counter = dim; e < dim; ++e)
        for (unsigned int d = e + 1; d < dim; ++d, ++counter)
          {
            const VectorizedArrayType value =
              sym_grad.access_raw_entry(counter) * JxW;
            gradients[e * nqp_d + d] = value * jac[d][d];
            gradients[d * nqp_d + e] = value * jac[e][e];
          }
    }
  // general/affine cell type
  else
    {
      const VectorizedArrayType JxW =
        this->cell_type == internal::MatrixFreeFunctions::general ?
          this->J_value[q_point] :
          this->J_value[0] * this->quadrature_weights[q_point];
      const Tensor<2, dim, VectorizedArrayType> jac =
        this->cell_type == internal::MatrixFreeFunctions::general ?
          this->jacobian[q_point] :
          this->jacobian[0];
      VectorizedArrayType weighted[dim][dim];
      for (unsigned int i = 0; i < dim; ++i)
        weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
      for (unsigned int i = 0, counter = dim; i < dim; ++i)
        for (unsigned int j = i + 1; j < dim; ++j, ++counter)
          {
            const VectorizedArrayType value =
              sym_grad.access_raw_entry(counter) * JxW;
            weighted[i][j] = value;
            weighted[j][i] = value;
          }
      for (unsigned int comp = 0; comp < dim; ++comp)
        for (unsigned int d = 0; d < dim; ++d)
          {
            VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
            for (unsigned int e = 1; e < dim; ++e)
              new_val += jac[e][d] * weighted[comp][e];
            gradients[comp * nqp_d + d] = new_val;
          }
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <int, typename>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_curl(
    const typename FEEvaluationBase<dim,
                                    n_components_,
                                    Number,
                                    is_face,
                                    VectorizedArrayType>::curl_type curl,
    const unsigned int                                              q_point)
{
  static_assert(n_components == dim,
                "Do not try to modify the default template parameters used for"
                " selectively enabling this function via std::enable_if!");
 
  Tensor<2, dim, VectorizedArrayType> grad;
 
  if (this->data->element_type ==
      internal::MatrixFreeFunctions::ElementType::tensor_nedelec)
    {
      curl_type curl_temp;
 
      // Piola transform is required
      AssertIndexRange(q_point, this->n_quadrature_points);
      Assert(this->J_value != nullptr,
             internal::ExcMatrixFreeAccessToUninitializedMappingField(
               "update_value"));
#  ifdef DEBUG
      Assert(this->is_reinitialized, ExcNotInitialized());
      this->gradients_quad_submitted = true;
#  endif
 
      const std::size_t    nqp_d     = this->n_quadrature_points * dim;
      VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
 
      if (!is_face &&
          this->cell_type == internal::MatrixFreeFunctions::cartesian)
        {
          const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
          const VectorizedArrayType weight = this->quadrature_weights[q_point];
 
          if constexpr (dim == 2)
            {
              curl_temp = curl * weight;
            }
          else if constexpr (dim == 3)
            {
              for (unsigned int comp = 0; comp < dim; ++comp)
                {
                  curl_temp[comp] = curl[comp] * jac[comp][comp] * weight;
                }
            }
          else
            DEAL_II_NOT_IMPLEMENTED();
        }
      else
        {
          // Affine or general cell
          const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
            (this->cell_type > internal::MatrixFreeFunctions::affine) ?
              this->jacobian[q_point] :
              this->jacobian[0];
 
          // Derivatives are reordered for faces. Need to take this into
          // account and 1/inv_det != J_value for faces
          const VectorizedArrayType fac =
            (!is_face) ?
              this->quadrature_weights[q_point] :
              (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
                  this->J_value[q_point] :
                  this->J_value[0] * this->quadrature_weights[q_point]) *
               ((dim == 2 && this->get_face_no() < 2) ?
                  -determinant(inv_t_jac) :
                  determinant(inv_t_jac)));
          const Tensor<2, dim, VectorizedArrayType> jac =
            (this->cell_type > internal::MatrixFreeFunctions::affine) ?
              transpose(invert(inv_t_jac)) :
              this->jacobian[1];
 
          if constexpr (dim == 2)
            {
              curl_temp = curl * fac;
            }
          else if constexpr (dim == 3)
            {
              // J^T * u * factor
              for (unsigned int comp = 0; comp < n_components; ++comp)
                {
                  curl_temp[comp] = curl[0] * jac[0][comp];
                  for (unsigned int e = 1; e < dim; ++e)
                    curl_temp[comp] += curl[e] * jac[e][comp];
                  curl_temp[comp] *= fac;
                }
            }
          else
            DEAL_II_NOT_IMPLEMENTED();
        }
      switch (dim)
        {
          case 2:
            grad[1][0] = curl_temp[0];
            grad[0][1] = -curl_temp[0];
            grad[0][0] = 0;
            grad[1][1] = 0;
            break;
          case 3:
            grad[2][1] = curl_temp[0];
            grad[1][2] = -curl_temp[0];
            grad[0][2] = curl_temp[1];
            grad[2][0] = -curl_temp[1];
            grad[1][0] = curl_temp[2];
            grad[0][1] = -curl_temp[2];
            grad[0][0] = 0;
            grad[1][1] = 0;
            grad[2][2] = 0;
            break;
          default:
            DEAL_II_NOT_IMPLEMENTED();
        }
      for (unsigned int d = 0; d < dim; ++d)
        {
          for (unsigned int comp = 0; comp < n_components; ++comp)
            gradients[comp * nqp_d + d] = grad[comp][d];
        }
      return;
    }
 
  if constexpr (dim == 2)
    {
      grad[1][0] = curl;
      grad[0][1] = -curl;
    }
  else if constexpr (dim == 3)
    {
      grad[2][1] = curl[0];
      grad[1][2] = -curl[0];
      grad[0][2] = curl[1];
      grad[2][0] = -curl[1];
      grad[1][0] = curl[2];
      grad[0][1] = -curl[2];
    }
  else
    DEAL_II_NOT_IMPLEMENTED();
 
  submit_gradient(grad, q_point);
}
 
 
 
/*-------------------------- FEEvaluation -----------------------------------*/
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline FEEvaluation<dim,
                    fe_degree,
                    n_q_points_1d,
                    n_components_,
                    Number,
                    VectorizedArrayType>::
  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
               const unsigned int                                  fe_no,
               const unsigned int quadrature_index,
               const unsigned int first_selected_component,
               const unsigned int active_fe_index,
               const unsigned int active_quad_index)
  : BaseClass(matrix_free,
              fe_no,
              first_selected_component,
              quadrature_index,
              fe_degree,
              static_n_q_points,
              true /*note: this is not a face*/,
              active_fe_index,
              active_quad_index,
              numbers::invalid_unsigned_int /*face_type*/)
  , dofs_per_component(this->data->dofs_per_component_on_cell)
  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
  , n_q_points(this->data->n_q_points)
{
  check_template_arguments(fe_no, 0);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline FEEvaluation<dim,
                    fe_degree,
                    n_q_points_1d,
                    n_components_,
                    Number,
                    VectorizedArrayType>::
  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
               const std::pair<unsigned int, unsigned int>        &range,
               const unsigned int dof_handler_index,
               const unsigned int quadrature_index,
               const unsigned int first_selected_component)
  : FEEvaluation(matrix_free,
                 dof_handler_index,
                 quadrature_index,
                 first_selected_component,
                 matrix_free.get_cell_active_fe_index(range, dof_handler_index))
{}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline FEEvaluation<dim,
                    fe_degree,
                    n_q_points_1d,
                    n_components_,
                    Number,
                    VectorizedArrayType>::
  FEEvaluation(const Mapping<dim>       &mapping,
               const FiniteElement<dim> &fe,
               const Quadrature<1>      &quadrature,
               const UpdateFlags         update_flags,
               const unsigned int        first_selected_component)
  : BaseClass(mapping,
              fe,
              quadrature,
              update_flags,
              first_selected_component,
              nullptr)
  , dofs_per_component(this->data->dofs_per_component_on_cell)
  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
  , n_q_points(this->data->n_q_points)
{
  check_template_arguments(numbers::invalid_unsigned_int, 0);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline FEEvaluation<dim,
                    fe_degree,
                    n_q_points_1d,
                    n_components_,
                    Number,
                    VectorizedArrayType>::
  FEEvaluation(const FiniteElement<dim> &fe,
               const Quadrature<1>      &quadrature,
               const UpdateFlags         update_flags,
               const unsigned int        first_selected_component)
  : BaseClass(StaticMappingQ1<dim>::mapping,
              fe,
              quadrature,
              update_flags,
              first_selected_component,
              nullptr)
  , dofs_per_component(this->data->dofs_per_component_on_cell)
  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
  , n_q_points(this->data->n_q_points)
{
  check_template_arguments(numbers::invalid_unsigned_int, 0);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline FEEvaluation<dim,
                    fe_degree,
                    n_q_points_1d,
                    n_components_,
                    Number,
                    VectorizedArrayType>::
  FEEvaluation(const FiniteElement<dim>                                &fe,
               const FEEvaluationData<dim, VectorizedArrayType, false> &other,
               const unsigned int first_selected_component)
  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
              fe,
              other.mapped_geometry->get_quadrature(),
              other.mapped_geometry->get_fe_values().get_update_flags(),
              first_selected_component,
              &other)
  , dofs_per_component(this->data->dofs_per_component_on_cell)
  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
  , n_q_points(this->data->n_q_points)
{
  check_template_arguments(numbers::invalid_unsigned_int, 0);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline FEEvaluation<dim,
                    fe_degree,
                    n_q_points_1d,
                    n_components_,
                    Number,
                    VectorizedArrayType>::FEEvaluation(const FEEvaluation
                                                         &other)
  : BaseClass(other)
  , dofs_per_component(this->data->dofs_per_component_on_cell)
  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
  , n_q_points(this->data->n_q_points)
{
  check_template_arguments(numbers::invalid_unsigned_int, 0);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline FEEvaluation<dim,
                    fe_degree,
                    n_q_points_1d,
                    n_components_,
                    Number,
                    VectorizedArrayType> &
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::operator=(const FEEvaluation &other)
{
  BaseClass::operator=(other);
  check_template_arguments(numbers::invalid_unsigned_int, 0);
  return *this;
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  check_template_arguments(const unsigned int dof_handler_index,
                           const unsigned int first_selected_component)
{
  (void)dof_handler_index;
  (void)first_selected_component;
 
  Assert(
    this->data->dofs_per_component_on_cell > 0,
    ExcMessage(
      "There is nothing useful you can do with an FEEvaluation object with "
      "FE_Nothing, i.e., without DoFs! If you have passed to "
      "MatrixFree::reinit() a collection of finite elements also containing "
      "FE_Nothing, please check - before creating FEEvaluation - the category "
      "of the current range by calling either "
      "MatrixFree::get_cell_range_category(range) or "
      "MatrixFree::get_face_range_category(range). The returned category "
      "is the index of the active FE, which you can use to exclude "
      "FE_Nothing."));
 
  if constexpr (running_in_debug_mode())
    {
      // print error message when the dimensions do not match. Propose a
      // possible fix
      if ((static_cast<unsigned int>(fe_degree) !=
             numbers::invalid_unsigned_int &&
           static_cast<unsigned int>(fe_degree) !=
             this->data->data.front().fe_degree) ||
          n_q_points != this->n_quadrature_points)
        {
          std::string message =
            "-------------------------------------------------------\n";
          message +=
            "Illegal arguments in constructor/wrong template arguments!\n";
          message += "    Called -->   FEEvaluation<dim,";
          message += Utilities::int_to_string(fe_degree) + ",";
          message += Utilities::int_to_string(n_q_points_1d);
          message += "," + Utilities::int_to_string(n_components);
          message += ",Number>(data";
          if (first_selected_component != numbers::invalid_unsigned_int)
            {
              message +=
                ", " + Utilities::int_to_string(dof_handler_index) + ", ";
              message +=
                Utilities::int_to_string(this->quadrature_index) + ", ";
              message += Utilities::int_to_string(first_selected_component);
            }
          message += ")\n";
 
          // check whether some other vector component has the correct number of
          // points
          unsigned int proposed_dof_comp  = numbers::invalid_unsigned_int,
                       proposed_fe_comp   = numbers::invalid_unsigned_int,
                       proposed_quad_comp = numbers::invalid_unsigned_int;
          if (dof_handler_index != numbers::invalid_unsigned_int)
            {
              if (static_cast<unsigned int>(fe_degree) ==
                  this->data->data.front().fe_degree)
                {
                  proposed_dof_comp = dof_handler_index;
                  proposed_fe_comp  = first_selected_component;
                }
              else
                for (unsigned int no = 0;
                     no < this->matrix_free->n_components();
                     ++no)
                  for (unsigned int nf = 0;
                       nf < this->matrix_free->n_base_elements(no);
                       ++nf)
                    if (this->matrix_free
                          ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
                          .data.front()
                          .fe_degree == static_cast<unsigned int>(fe_degree))
                      {
                        proposed_dof_comp = no;
                        proposed_fe_comp  = nf;
                        break;
                      }
              if (n_q_points ==
                  this->mapping_data->descriptor[this->active_quad_index]
                    .n_q_points)
                proposed_quad_comp = this->quadrature_index;
              else
                for (unsigned int no = 0;
                     no <
                     this->matrix_free->get_mapping_info().cell_data.size();
                     ++no)
                  if (this->matrix_free->get_mapping_info()
                        .cell_data[no]
                        .descriptor[this->active_quad_index]
                        .n_q_points == n_q_points)
                    {
                      proposed_quad_comp = no;
                      break;
                    }
            }
          if (proposed_dof_comp != numbers::invalid_unsigned_int &&
              proposed_quad_comp != numbers::invalid_unsigned_int)
            {
              if (proposed_dof_comp != first_selected_component)
                message += "Wrong vector component selection:\n";
              else
                message += "Wrong quadrature formula selection:\n";
              message += "    Did you mean FEEvaluation<dim,";
              message += Utilities::int_to_string(fe_degree) + ",";
              message += Utilities::int_to_string(n_q_points_1d);
              message += "," + Utilities::int_to_string(n_components);
              message += ",Number>(data";
              if (dof_handler_index != numbers::invalid_unsigned_int)
                {
                  message +=
                    ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
                  message +=
                    Utilities::int_to_string(proposed_quad_comp) + ", ";
                  message += Utilities::int_to_string(proposed_fe_comp);
                }
              message += ")?\n";
              std::string correct_pos;
              if (proposed_dof_comp != dof_handler_index)
                correct_pos = " ^ ";
              else
                correct_pos = "   ";
              if (proposed_quad_comp != this->quadrature_index)
                correct_pos += " ^ ";
              else
                correct_pos += "   ";
              if (proposed_fe_comp != first_selected_component)
                correct_pos += " ^\n";
              else
                correct_pos += "  \n";
              message +=
                "                                                     " +
                correct_pos;
            }
          // ok, did not find the numbers specified by the template arguments in
          // the given list. Suggest correct template arguments
          const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
            std::pow(1.001 * this->n_quadrature_points, 1. / dim));
          message += "Wrong template arguments:\n";
          message += "    Did you mean FEEvaluation<dim,";
          message +=
            Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
          message += Utilities::int_to_string(proposed_n_q_points_1d);
          message += "," + Utilities::int_to_string(n_components);
          message += ",Number>(data";
          if (dof_handler_index != numbers::invalid_unsigned_int)
            {
              message +=
                ", " + Utilities::int_to_string(dof_handler_index) + ", ";
              message += Utilities::int_to_string(this->quadrature_index);
              message +=
                ", " + Utilities::int_to_string(first_selected_component);
            }
          message += ")?\n";
          std::string correct_pos;
          if (this->data->data.front().fe_degree !=
              static_cast<unsigned int>(fe_degree))
            correct_pos = " ^";
          else
            correct_pos = "  ";
          if (proposed_n_q_points_1d != n_q_points_1d)
            correct_pos += " ^\n";
          else
            correct_pos += "  \n";
          message += "                                 " + correct_pos;
 
          Assert(static_cast<unsigned int>(fe_degree) ==
                     this->data->data.front().fe_degree &&
                   n_q_points == this->n_quadrature_points,
                 ExcMessage(message));
        }
      if (dof_handler_index != numbers::invalid_unsigned_int)
        AssertDimension(
          n_q_points,
          this->mapping_data->descriptor[this->active_quad_index].n_q_points);
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::reinit(const unsigned int cell_index)
{
  Assert(this->matrix_free != nullptr,
         ExcMessage("FEEvaluation was initialized without a matrix-free object."
                    " Integer indexing is not possible."));
 
  Assert(this->dof_info != nullptr, ExcNotInitialized());
  Assert(this->mapping_data != nullptr, ExcNotInitialized());
  this->cell = cell_index;
  this->cell_type =
    this->matrix_free->get_mapping_info().get_cell_type(cell_index);
 
  const unsigned int offsets =
    this->mapping_data->data_index_offsets[cell_index];
  this->jacobian = &this->mapping_data->jacobians[0][offsets];
  this->J_value  = &this->mapping_data->JxW_values[offsets];
  if (!this->mapping_data->jacobian_gradients[0].empty())
    {
      this->jacobian_gradients =
        this->mapping_data->jacobian_gradients[0].data() + offsets;
      this->jacobian_gradients_non_inverse =
        this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
    }
 
  if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) == n_lanes)
    {
      DEAL_II_OPENMP_SIMD_PRAGMA
      for (unsigned int i = 0; i < n_lanes; ++i)
        this->cell_ids[i] = cell_index * n_lanes + i;
    }
  else
    {
      unsigned int i = 0;
      for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
           ++i)
        this->cell_ids[i] = cell_index * n_lanes + i;
      for (; i < n_lanes; ++i)
        this->cell_ids[i] = numbers::invalid_unsigned_int;
    }
 
  if (this->mapping_data->quadrature_points.empty() == false)
    this->quadrature_points =
      &this->mapping_data->quadrature_points
         [this->mapping_data->quadrature_point_offsets[this->cell]];
 
  if constexpr (running_in_debug_mode())
    {
      this->is_reinitialized           = true;
      this->dof_values_initialized     = false;
      this->values_quad_initialized    = false;
      this->gradients_quad_initialized = false;
      this->hessians_quad_initialized  = false;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::reinit(const std::array<unsigned int,
                                                           n_lanes> &cell_ids)
{
  Assert(this->dof_info != nullptr, ExcNotInitialized());
  Assert(this->mapping_data != nullptr, ExcNotInitialized());
 
  this->cell     = numbers::invalid_unsigned_int;
  this->cell_ids = cell_ids;
 
  // determine type of cell batch
  this->cell_type = internal::MatrixFreeFunctions::GeometryType::cartesian;
 
  for (unsigned int v = 0; v < n_lanes; ++v)
    {
      const unsigned int cell_index = cell_ids[v];
 
      if (cell_index == numbers::invalid_unsigned_int)
        continue;
 
      this->cell_type =
        std::max(this->cell_type,
                 this->matrix_free->get_mapping_info().get_cell_type(
                   cell_index / n_lanes));
    }
 
  // allocate memory for internal data storage
  if (this->mapped_geometry == nullptr)
    this->mapped_geometry =
      std::make_shared<internal::MatrixFreeFunctions::
                         MappingDataOnTheFly<dim, VectorizedArrayType>>();
 
  auto &mapping_storage = this->mapped_geometry->get_data_storage();
 
  auto &this_jacobian_data           = mapping_storage.jacobians[0];
  auto &this_J_value_data            = mapping_storage.JxW_values;
  auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
  auto &this_jacobian_gradients_non_inverse_data =
    mapping_storage.jacobian_gradients_non_inverse[0];
  auto &this_quadrature_points_data = mapping_storage.quadrature_points;
 
  if (this->cell_type <= internal::MatrixFreeFunctions::GeometryType::affine)
    {
      if (this_jacobian_data.size() != 2)
        this_jacobian_data.resize_fast(2);
 
      if (this_J_value_data.size() != 1)
        this_J_value_data.resize_fast(1);
 
      const auto &update_flags_cells =
        this->matrix_free->get_mapping_info().update_flags_cells;
 
      if (update_flags_cells & update_jacobian_grads &&
          this_jacobian_gradients_data.size() != 1)
        {
          this_jacobian_gradients_data.resize_fast(1);
          this_jacobian_gradients_non_inverse_data.resize_fast(1);
        }
 
      if (update_flags_cells & update_quadrature_points &&
          this_quadrature_points_data.size() != 1)
        this_quadrature_points_data.resize_fast(1);
    }
  else
    {
      if (this_jacobian_data.size() != this->n_quadrature_points)
        this_jacobian_data.resize_fast(this->n_quadrature_points);
 
      if (this_J_value_data.size() != this->n_quadrature_points)
        this_J_value_data.resize_fast(this->n_quadrature_points);
 
      const auto &update_flags_cells =
        this->matrix_free->get_mapping_info().update_flags_cells;
 
      if (update_flags_cells & update_jacobian_grads &&
          this_jacobian_gradients_data.size() != this->n_quadrature_points)
        {
          this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
          this_jacobian_gradients_non_inverse_data.resize_fast(
            this->n_quadrature_points);
        }
 
      if (update_flags_cells & update_quadrature_points &&
          this_quadrature_points_data.size() != this->n_quadrature_points)
        this_quadrature_points_data.resize_fast(this->n_quadrature_points);
    }
 
  // set pointers to internal data storage
  this->jacobian           = this_jacobian_data.data();
  this->J_value            = this_J_value_data.data();
  this->jacobian_gradients = this_jacobian_gradients_data.data();
  this->jacobian_gradients_non_inverse =
    this_jacobian_gradients_non_inverse_data.data();
  this->quadrature_points = this_quadrature_points_data.data();
 
  // fill internal data storage lane by lane
  for (unsigned int v = 0; v < n_lanes; ++v)
    {
      const unsigned int cell_index = cell_ids[v];
 
      if (cell_index == numbers::invalid_unsigned_int)
        continue;
 
      const unsigned int cell_batch_index = cell_index / n_lanes;
      const unsigned int offsets =
        this->mapping_data->data_index_offsets[cell_batch_index];
      const unsigned int lane = cell_index % n_lanes;
 
      if (this->cell_type <=
          internal::MatrixFreeFunctions::GeometryType::affine)
        {
          // case that all cells are Cartesian or affine
          for (unsigned int q = 0; q < 2; ++q)
            for (unsigned int i = 0; i < dim; ++i)
              for (unsigned int j = 0; j < dim; ++j)
                this_jacobian_data[q][i][j][v] =
                  this->mapping_data->jacobians[0][offsets + q][i][j][lane];
 
          const unsigned int q = 0;
 
          this_J_value_data[q][v] =
            this->mapping_data->JxW_values[offsets + q][lane];
 
          const auto &update_flags_cells =
            this->matrix_free->get_mapping_info().update_flags_cells;
 
          if (update_flags_cells & update_jacobian_grads)
            {
              for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
                for (unsigned int j = 0; j < dim; ++j)
                  this_jacobian_gradients_data[q][i][j][v] =
                    this->mapping_data
                      ->jacobian_gradients[0][offsets + q][i][j][lane];
 
              for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
                for (unsigned int j = 0; j < dim; ++j)
                  this_jacobian_gradients_non_inverse_data[q][i][j][v] =
                    this->mapping_data
                      ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
                                                      [lane];
            }
 
          if (update_flags_cells & update_quadrature_points)
            for (unsigned int i = 0; i < dim; ++i)
              this_quadrature_points_data[q][i][v] =
                this->mapping_data->quadrature_points
                  [this->mapping_data
                     ->quadrature_point_offsets[cell_batch_index] +
                   q][i][lane];
        }
      else
        {
          // general case that at least one cell is not Cartesian or affine
          const auto cell_type =
            this->matrix_free->get_mapping_info().get_cell_type(
              cell_batch_index);
 
          for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
            {
              const unsigned int q_src =
                (cell_type <=
                 internal::MatrixFreeFunctions::GeometryType::affine) ?
                  0 :
                  q;
 
              this_J_value_data[q][v] =
                this->mapping_data->JxW_values[offsets + q_src][lane];
 
              for (unsigned int i = 0; i < dim; ++i)
                for (unsigned int j = 0; j < dim; ++j)
                  this_jacobian_data[q][i][j][v] =
                    this->mapping_data
                      ->jacobians[0][offsets + q_src][i][j][lane];
 
              const auto &update_flags_cells =
                this->matrix_free->get_mapping_info().update_flags_cells;
 
              if (update_flags_cells & update_jacobian_grads)
                {
                  for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
                    for (unsigned int j = 0; j < dim; ++j)
                      this_jacobian_gradients_data[q][i][j][v] =
                        this->mapping_data
                          ->jacobian_gradients[0][offsets + q_src][i][j][lane];
 
                  for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
                    for (unsigned int j = 0; j < dim; ++j)
                      this_jacobian_gradients_non_inverse_data[q][i][j][v] =
                        this->mapping_data
                          ->jacobian_gradients_non_inverse[0][offsets + q_src]
                                                          [i][j][lane];
                }
 
              if (update_flags_cells & update_quadrature_points)
                {
                  if (cell_type <=
                      internal::MatrixFreeFunctions::GeometryType::affine)
                    {
                      // affine case: quadrature points are not available but
                      // have to be computed from the corner point and the
                      // Jacobian
                      Point<dim, VectorizedArrayType> point =
                        this->mapping_data->quadrature_points
                          [this->mapping_data
                             ->quadrature_point_offsets[cell_batch_index] +
                           0];
 
                      const Tensor<2, dim, VectorizedArrayType> &jac =
                        this->mapping_data->jacobians[0][offsets + 1];
                      if (cell_type == internal::MatrixFreeFunctions::cartesian)
                        for (unsigned int d = 0; d < dim; ++d)
                          point[d] +=
                            jac[d][d] *
                            static_cast<Number>(
                              this->descriptor->quadrature.point(q)[d]);
                      else
                        for (unsigned int d = 0; d < dim; ++d)
                          for (unsigned int e = 0; e < dim; ++e)
                            point[d] +=
                              jac[d][e] *
                              static_cast<Number>(
                                this->descriptor->quadrature.point(q)[e]);
 
                      for (unsigned int i = 0; i < dim; ++i)
                        this_quadrature_points_data[q][i][v] = point[i][lane];
                    }
                  else
                    {
                      // general case: quadrature points are available
                      for (unsigned int i = 0; i < dim; ++i)
                        this_quadrature_points_data[q][i][v] =
                          this->mapping_data->quadrature_points
                            [this->mapping_data
                               ->quadrature_point_offsets[cell_batch_index] +
                             q][i][lane];
                    }
                }
            }
        }
    }
 
  if constexpr (running_in_debug_mode())
    {
      this->is_reinitialized           = true;
      this->dof_values_initialized     = false;
      this->values_quad_initialized    = false;
      this->gradients_quad_initialized = false;
      this->hessians_quad_initialized  = false;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
template <bool level_dof_access>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  reinit(const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell)
{
  Assert(this->matrix_free == nullptr,
         ExcMessage("Cannot use initialization from cell iterator if "
                    "initialized from MatrixFree object. Use variant for "
                    "on the fly computation with arguments as for FEValues "
                    "instead"));
  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
  this->mapped_geometry->reinit(
    static_cast<typename Triangulation<dim>::cell_iterator>(cell));
  this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
  if (level_dof_access)
    cell->get_mg_dof_indices(this->local_dof_indices);
  else
    cell->get_dof_indices(this->local_dof_indices);
 
  if constexpr (running_in_debug_mode())
    {
      this->is_reinitialized = true;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  reinit(const typename Triangulation<dim>::cell_iterator &cell)
{
  Assert(this->matrix_free == 0,
         ExcMessage("Cannot use initialization from cell iterator if "
                    "initialized from MatrixFree object. Use variant for "
                    "on the fly computation with arguments as for FEValues "
                    "instead"));
  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
  this->mapped_geometry->reinit(cell);
 
  if constexpr (running_in_debug_mode())
    {
      this->is_reinitialized = true;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->dof_values_initialized == true,
             internal::ExcAccessToUninitializedField());
    }
  evaluate(this->values_dofs, evaluation_flags);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  evaluate(const VectorizedArrayType             *values_array,
           const EvaluationFlags::EvaluationFlags evaluation_flag)
{
  const bool hessians_on_general_cells =
    evaluation_flag & EvaluationFlags::hessians &&
    (this->cell_type > internal::MatrixFreeFunctions::affine);
  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
  if (hessians_on_general_cells)
    evaluation_flag_actual |= EvaluationFlags::gradients;
 
  if (this->data->element_type ==
        internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
      evaluation_flag & EvaluationFlags::gradients &&
      (this->cell_type > internal::MatrixFreeFunctions::affine))
    evaluation_flag_actual |= EvaluationFlags::values;
 
  if constexpr (fe_degree > -1)
    {
      internal::FEEvaluationImplSelector<dim, VectorizedArrayType, false>::
        template run<fe_degree, n_q_points_1d>(n_components,
                                               evaluation_flag_actual,
                                               values_array,
                                               *this);
    }
  else
    {
      internal::FEEvaluationFactory<dim, VectorizedArrayType>::evaluate(
        n_components,
        evaluation_flag_actual,
        const_cast<VectorizedArrayType *>(values_array),
        *this);
    }
 
  if constexpr (running_in_debug_mode())
    {
      this->values_quad_initialized =
        evaluation_flag_actual & EvaluationFlags::values;
      this->gradients_quad_initialized =
        evaluation_flag_actual & EvaluationFlags::gradients;
      this->hessians_quad_initialized =
        evaluation_flag_actual & EvaluationFlags::hessians;
    }
}
 
 
namespace internal
{
  template <typename Number,
            typename VectorizedArrayType,
            typename VectorType,
            typename EvaluatorType,
            std::enable_if_t<internal::has_begin<VectorType> &&
                               !IsBlockVector<VectorType>::value,
                             VectorType> * = nullptr>
  VectorizedArrayType *
  check_vector_access_inplace(const EvaluatorType &fe_eval, VectorType &vector)
  {
    // for user-defined cell batches this functionality is not supported
    if (fe_eval.get_current_cell_index() == numbers::invalid_unsigned_int)
      return nullptr;
 
    const unsigned int cell     = fe_eval.get_cell_or_face_batch_id();
    const auto        &dof_info = fe_eval.get_dof_info();
 
    // If the index storage is interleaved and contiguous and the vector
    // storage has the correct alignment, we can directly pass the pointer
    // into the vector to the evaluate() and integrate() calls, without
    // reading the vector entries into a separate data field. This saves some
    // operations.
    if (std::is_same_v<typename VectorType::value_type, Number> &&
        dof_info.index_storage_variants
            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] ==
          internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
            interleaved_contiguous &&
        reinterpret_cast<std::size_t>(
          vector.begin() +
          dof_info.dof_indices_contiguous
            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
            [cell * VectorizedArrayType::size()]) %
            sizeof(VectorizedArrayType) ==
          0)
      {
        return reinterpret_cast<VectorizedArrayType *>(
          vector.begin() +
          dof_info.dof_indices_contiguous
            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
            [cell * VectorizedArrayType::size()] +
          dof_info.component_dof_indices_offset
              [fe_eval.get_active_fe_index()]
              [fe_eval.get_first_selected_component()] *
            VectorizedArrayType::size());
      }
    else
      return nullptr;
  }
 
  template <typename Number,
            typename VectorizedArrayType,
            typename VectorType,
            typename EvaluatorType,
            std::enable_if_t<!internal::has_begin<VectorType> ||
                               IsBlockVector<VectorType>::value,
                             VectorType> * = nullptr>
  VectorizedArrayType *
  check_vector_access_inplace(const EvaluatorType &, VectorType &)
  {
    return nullptr;
  }
} // namespace internal
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  gather_evaluate(const VectorType                      &input_vector,
                  const EvaluationFlags::EvaluationFlags evaluation_flag)
{
  const VectorizedArrayType *src_ptr =
    internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
      *this, input_vector);
  if (src_ptr != nullptr)
    evaluate(src_ptr, evaluation_flag);
  else
    {
      this->read_dof_values(input_vector);
      evaluate(this->begin_dof_values(), evaluation_flag);
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  integrate(const EvaluationFlags::EvaluationFlags integration_flag)
{
  integrate(integration_flag, this->values_dofs);
 
  if constexpr (running_in_debug_mode())
    {
      this->dof_values_initialized = true;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
            VectorizedArrayType                   *values_array,
            const bool                             sum_into_values_array)
{
  if constexpr (running_in_debug_mode())
    {
      if (integration_flag & EvaluationFlags::values)
        Assert(this->values_quad_submitted == true,
               internal::ExcAccessToUninitializedField());
      if (integration_flag & EvaluationFlags::gradients)
        Assert(this->gradients_quad_submitted == true,
               internal::ExcAccessToUninitializedField());
      if ((integration_flag & EvaluationFlags::hessians) != 0u)
        Assert(this->hessians_quad_submitted == true,
               internal::ExcAccessToUninitializedField());
    }
  Assert(this->matrix_free != nullptr ||
           this->mapped_geometry->is_initialized(),
         ExcNotInitialized());
 
  Assert(
    (integration_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
                          EvaluationFlags::hessians)) == 0,
    ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, and "
               "EvaluationFlags::hessians are supported."));
 
  EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
  if (integration_flag & EvaluationFlags::hessians &&
      (this->cell_type > internal::MatrixFreeFunctions::affine))
    {
      unsigned int size = n_components * dim * n_q_points;
      if ((integration_flag & EvaluationFlags::gradients) != 0u)
        {
          for (unsigned int i = 0; i < size; ++i)
            this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
        }
      else
        {
          for (unsigned int i = 0; i < size; ++i)
            this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
          integration_flag_actual |= EvaluationFlags::gradients;
        }
    }
 
  if (n_components == dim &&
      this->data->element_type ==
        internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
      integration_flag & EvaluationFlags::gradients &&
      this->cell_type > internal::MatrixFreeFunctions::affine &&
      this->divergence_is_requested == false)
    {
      unsigned int size = n_components * n_q_points;
      if ((integration_flag & EvaluationFlags::values) != 0u)
        {
          for (unsigned int i = 0; i < size; ++i)
            this->values_quad[i] += this->values_from_gradients_quad[i];
        }
      else
        {
          for (unsigned int i = 0; i < size; ++i)
            this->values_quad[i] = this->values_from_gradients_quad[i];
          integration_flag_actual |= EvaluationFlags::values;
        }
    }
 
  if constexpr (fe_degree > -1)
    {
      internal::FEEvaluationImplSelector<dim, VectorizedArrayType, true>::
        template run<fe_degree, n_q_points_1d>(n_components,
                                               integration_flag_actual,
                                               values_array,
                                               *this,
                                               sum_into_values_array);
    }
  else
    {
      internal::FEEvaluationFactory<dim, VectorizedArrayType>::integrate(
        n_components,
        integration_flag_actual,
        values_array,
        *this,
        sum_into_values_array);
    }
 
  if constexpr (running_in_debug_mode())
    {
      this->dof_values_initialized = true;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
                    VectorType                            &destination)
{
  VectorizedArrayType *dst_ptr =
    internal::check_vector_access_inplace<Number, VectorizedArrayType>(
      *this, destination);
  if (dst_ptr != nullptr)
    integrate(integration_flag, dst_ptr, true);
  else
    {
      integrate(integration_flag, this->begin_dof_values());
      this->distribute_local_to_global(destination);
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::dof_indices() const
{
  return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
    0U, dofs_per_cell);
}
 
 
 
/*-------------------------- FEFaceEvaluation ---------------------------*/
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline FEFaceEvaluation<dim,
                        fe_degree,
                        n_q_points_1d,
                        n_components_,
                        Number,
                        VectorizedArrayType>::
  FEFaceEvaluation(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const bool                                          is_interior_face,
    const unsigned int                                  dof_handler_index,
    const unsigned int                                  quadrature_index,
    const unsigned int first_selected_component,
    const unsigned int active_fe_index,
    const unsigned int active_quad_index,
    const unsigned int face_type)
  : BaseClass(matrix_free,
              dof_handler_index,
              first_selected_component,
              quadrature_index,
              fe_degree,
              static_n_q_points,
              is_interior_face,
              active_fe_index,
              active_quad_index,
              face_type)
  , dofs_per_component(this->data->dofs_per_component_on_cell)
  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
  , n_q_points(this->n_quadrature_points)
{}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline FEFaceEvaluation<dim,
                        fe_degree,
                        n_q_points_1d,
                        n_components_,
                        Number,
                        VectorizedArrayType>::
  FEFaceEvaluation(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const std::pair<unsigned int, unsigned int>        &range,
    const bool                                          is_interior_face,
    const unsigned int                                  dof_handler_index,
    const unsigned int                                  quadrature_index,
    const unsigned int first_selected_component)
  : FEFaceEvaluation(matrix_free,
                     is_interior_face,
                     dof_handler_index,
                     quadrature_index,
                     first_selected_component,
                     matrix_free.get_face_active_fe_index(range,
                                                          is_interior_face),
                     numbers::invalid_unsigned_int,
                     matrix_free.get_face_info(range.first).face_type)
{}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::reinit(const unsigned int face_index)
{
  Assert(this->mapped_geometry == nullptr,
         ExcMessage("FEEvaluation was initialized without a matrix-free object."
                    " Integer indexing is not possible"));
  if (this->mapped_geometry != nullptr)
    return;
 
  this->cell = face_index;
  this->dof_access_index =
    this->is_interior_face() ?
      internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
      internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior;
  Assert(this->mapping_data != nullptr, ExcNotInitialized());
 
  if (face_index >=
        this->matrix_free->get_task_info().face_partition_data.back() &&
      face_index <
        this->matrix_free->get_task_info().boundary_partition_data.back())
    Assert(this->is_interior_face(),
           ExcMessage(
             "Boundary faces do not have a neighbor. When looping over "
             "boundary faces use FEFaceEvaluation with the parameter "
             "is_interior_face set to true. "));
 
  this->reinit_face(this->matrix_free->get_face_info(face_index));
 
  unsigned int i = 0;
  for (; i < this->matrix_free->n_active_entries_per_face_batch(this->cell);
       ++i)
    this->face_ids[i] = face_index * n_lanes + i;
  for (; i < n_lanes; ++i)
    this->face_ids[i] = numbers::invalid_unsigned_int;
 
  this->cell_type = this->matrix_free->get_mapping_info().face_type[face_index];
  const unsigned int offsets =
    this->mapping_data->data_index_offsets[face_index];
  this->J_value        = &this->mapping_data->JxW_values[offsets];
  this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
  this->jacobian =
    &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
  this->normal_x_jacobian =
    &this->mapping_data
       ->normals_times_jacobians[!this->is_interior_face()][offsets];
  this->jacobian_gradients =
    this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
    offsets;
  this->jacobian_gradients_non_inverse =
    this->mapping_data
      ->jacobian_gradients_non_inverse[!this->is_interior_face()]
      .data() +
    offsets;
 
  if (this->mapping_data->quadrature_point_offsets.empty() == false)
    {
      AssertIndexRange(this->cell,
                       this->mapping_data->quadrature_point_offsets.size());
      this->quadrature_points =
        this->mapping_data->quadrature_points.data() +
        this->mapping_data->quadrature_point_offsets[this->cell];
    }
 
  if constexpr (running_in_debug_mode())
    {
      this->is_reinitialized           = true;
      this->dof_values_initialized     = false;
      this->values_quad_initialized    = false;
      this->gradients_quad_initialized = false;
      this->hessians_quad_initialized  = false;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::reinit(const unsigned int cell_index,
                                              const unsigned int face_number)
{
  Assert(
    this->quadrature_index <
      this->matrix_free->get_mapping_info().face_data_by_cells.size(),
    ExcMessage(
      "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
  AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
  AssertIndexRange(cell_index,
                   this->matrix_free->get_mapping_info().cell_type.size());
  Assert(this->mapped_geometry == nullptr,
         ExcMessage("FEEvaluation was initialized without a matrix-free object."
                    " Integer indexing is not possible"));
  if (this->mapped_geometry != nullptr)
    return;
  Assert(this->matrix_free != nullptr, ExcNotInitialized());
 
  this->cell_type = this->matrix_free->get_mapping_info()
                      .faces_by_cells_type[cell_index][face_number];
  this->cell          = cell_index;
  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
  this->dof_access_index =
    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
 
  if (this->is_interior_face() == false)
    {
      // for this case, we need to look into the FaceInfo field that collects
      // information from both sides of a face once for the global mesh, and
      // pick the face id that is not the local one (cell_this).
      for (unsigned int i = 0; i < n_lanes; ++i)
        {
          // compute actual (non vectorized) cell ID
          const unsigned int cell_this = cell_index * n_lanes + i;
          // compute face ID
          unsigned int face_index =
            this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
                                                                  face_number,
                                                                  i);
 
          this->face_ids[i] = face_index;
 
          if (face_index == numbers::invalid_unsigned_int)
            {
              this->cell_ids[i]     = numbers::invalid_unsigned_int;
              this->face_numbers[i] = static_cast<std::uint8_t>(-1);
              this->face_orientations[i] =
                numbers::invalid_geometric_orientation;
              continue; // invalid face ID: no neighbor on boundary
            }
 
          const auto &faces =
            this->matrix_free->get_face_info(face_index / n_lanes);
          // get cell ID on both sides of face
          auto cell_m = faces.cells_interior[face_index % n_lanes];
          auto cell_p = faces.cells_exterior[face_index % n_lanes];
 
          const bool face_identifies_as_interior = cell_m != cell_this;
 
          Assert(cell_m == cell_this || cell_p == cell_this,
                 ExcInternalError());
 
          // compare the IDs with the given cell ID
          if (face_identifies_as_interior)
            {
              this->cell_ids[i]     = cell_m; // neighbor has the other ID
              this->face_numbers[i] = faces.interior_face_no;
            }
          else
            {
              this->cell_ids[i]     = cell_p;
              this->face_numbers[i] = faces.exterior_face_no;
            }
 
          const bool orientation_interior_face = faces.face_orientation >= 8;
          types::geometric_orientation face_orientation =
            faces.face_orientation % 8;
          if (face_identifies_as_interior != orientation_interior_face)
            {
              Assert(this->matrix_free->get_cell_iterator(cell_index, i)
                         ->reference_cell() ==
                       ReferenceCells::get_hypercube<dim>(),
                     ExcNotImplemented());
              face_orientation =
                ReferenceCells::get_hypercube<dim - 1>()
                  .get_inverse_combined_orientation(face_orientation);
            }
          this->face_orientations[i] = face_orientation;
        }
    }
  else
    {
      this->face_orientations[0] = numbers::default_geometric_orientation;
      this->face_numbers[0]      = face_number;
      if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) ==
          n_lanes)
        {
          DEAL_II_OPENMP_SIMD_PRAGMA
          for (unsigned int i = 0; i < n_lanes; ++i)
            this->cell_ids[i] = cell_index * n_lanes + i;
        }
      else
        {
          unsigned int i = 0;
          for (; i <
                 this->matrix_free->n_active_entries_per_cell_batch(this->cell);
               ++i)
            this->cell_ids[i] = cell_index * n_lanes + i;
          for (; i < n_lanes; ++i)
            this->cell_ids[i] = numbers::invalid_unsigned_int;
        }
      for (unsigned int i = 0; i < n_lanes; ++i)
        this->face_ids[i] =
          this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
                                                                face_number,
                                                                i);
    }
 
  const unsigned int offsets =
    this->matrix_free->get_mapping_info()
      .face_data_by_cells[this->quadrature_index]
      .data_index_offsets[cell_index * GeometryInfo<dim>::faces_per_cell +
                          face_number];
  AssertIndexRange(offsets,
                   this->matrix_free->get_mapping_info()
                     .face_data_by_cells[this->quadrature_index]
                     .JxW_values.size());
  this->J_value = &this->matrix_free->get_mapping_info()
                     .face_data_by_cells[this->quadrature_index]
                     .JxW_values[offsets];
  this->normal_vectors = &this->matrix_free->get_mapping_info()
                            .face_data_by_cells[this->quadrature_index]
                            .normal_vectors[offsets];
  this->jacobian = &this->matrix_free->get_mapping_info()
                      .face_data_by_cells[this->quadrature_index]
                      .jacobians[!this->is_interior_face()][offsets];
  this->normal_x_jacobian =
    &this->matrix_free->get_mapping_info()
       .face_data_by_cells[this->quadrature_index]
       .normals_times_jacobians[!this->is_interior_face()][offsets];
  this->jacobian_gradients =
    this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
    offsets;
  this->jacobian_gradients_non_inverse =
    this->mapping_data
      ->jacobian_gradients_non_inverse[!this->is_interior_face()]
      .data() +
    offsets;
 
  if (this->matrix_free->get_mapping_info()
        .face_data_by_cells[this->quadrature_index]
        .quadrature_point_offsets.empty() == false)
    {
      const unsigned int index =
        this->cell * GeometryInfo<dim>::faces_per_cell + this->face_numbers[0];
      AssertIndexRange(index,
                       this->matrix_free->get_mapping_info()
                         .face_data_by_cells[this->quadrature_index]
                         .quadrature_point_offsets.size());
      this->quadrature_points = this->matrix_free->get_mapping_info()
                                  .face_data_by_cells[this->quadrature_index]
                                  .quadrature_points.data() +
                                this->matrix_free->get_mapping_info()
                                  .face_data_by_cells[this->quadrature_index]
                                  .quadrature_point_offsets[index];
    }
 
  if constexpr (running_in_debug_mode())
    {
      this->is_reinitialized           = true;
      this->dof_values_initialized     = false;
      this->values_quad_initialized    = false;
      this->gradients_quad_initialized = false;
      this->hessians_quad_initialized  = false;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->dof_values_initialized, ExcNotInitialized());
    }
 
  evaluate(this->values_dofs, evaluation_flag);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  evaluate(const VectorizedArrayType             *values_array,
           const EvaluationFlags::EvaluationFlags evaluation_flag)
{
  Assert((evaluation_flag &
          ~(EvaluationFlags::values | EvaluationFlags::gradients |
            EvaluationFlags::hessians)) == 0,
         ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
                    "and EvaluationFlags::hessians are supported."));
 
  const bool hessians_on_general_cells =
    evaluation_flag & EvaluationFlags::hessians &&
    (this->cell_type > internal::MatrixFreeFunctions::affine);
  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
  if (hessians_on_general_cells)
    evaluation_flag_actual |= EvaluationFlags::gradients;
 
  if (this->data->element_type ==
        internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
      evaluation_flag & EvaluationFlags::gradients &&
      (this->cell_type > internal::MatrixFreeFunctions::affine))
    evaluation_flag_actual |= EvaluationFlags::values;
 
  if constexpr (fe_degree > -1)
    internal::FEFaceEvaluationImplEvaluateSelector<dim, VectorizedArrayType>::
      template run<fe_degree, n_q_points_1d>(n_components,
                                             evaluation_flag_actual,
                                             values_array,
                                             *this);
  else
    internal::FEFaceEvaluationFactory<dim, VectorizedArrayType>::evaluate(
      n_components, evaluation_flag_actual, values_array, *this);
 
  if constexpr (running_in_debug_mode())
    {
      this->values_quad_initialized =
        evaluation_flag_actual & EvaluationFlags::values;
      this->gradients_quad_initialized =
        evaluation_flag_actual & EvaluationFlags::gradients;
      this->hessians_quad_initialized =
        evaluation_flag_actual & EvaluationFlags::hessians;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  project_to_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
{
  if constexpr (running_in_debug_mode())
    {
      Assert(this->dof_values_initialized, ExcNotInitialized());
    }
 
  project_to_face(this->values_dofs, evaluation_flag);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  project_to_face(const VectorizedArrayType             *values_array,
                  const EvaluationFlags::EvaluationFlags evaluation_flag)
{
  Assert((evaluation_flag &
          ~(EvaluationFlags::values | EvaluationFlags::gradients |
            EvaluationFlags::hessians)) == 0,
         ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
                    "and EvaluationFlags::hessians are supported."));
 
  const bool hessians_on_general_cells =
    evaluation_flag & EvaluationFlags::hessians &&
    (this->cell_type > internal::MatrixFreeFunctions::affine);
  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
  if (hessians_on_general_cells)
    evaluation_flag_actual |= EvaluationFlags::gradients;
 
  if (this->data->element_type ==
        internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
      evaluation_flag & EvaluationFlags::gradients &&
      (this->cell_type > internal::MatrixFreeFunctions::affine))
    evaluation_flag_actual |= EvaluationFlags::values;
 
  if constexpr (fe_degree > -1)
    internal::FEFaceEvaluationImplProjectToFaceSelector<
      dim,
      VectorizedArrayType>::template run<fe_degree>(n_components,
                                                    evaluation_flag_actual,
                                                    values_array,
                                                    *this);
  else
    internal::FEFaceEvaluationFactory<dim, VectorizedArrayType>::
      project_to_face(n_components,
                      evaluation_flag_actual,
                      values_array,
                      *this);
 
  // face dofs initialized
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  evaluate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
{
  Assert((evaluation_flag &
          ~(EvaluationFlags::values | EvaluationFlags::gradients |
            EvaluationFlags::hessians)) == 0,
         ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
                    "and EvaluationFlags::hessians are supported."));
 
  const bool hessians_on_general_cells =
    evaluation_flag & EvaluationFlags::hessians &&
    (this->cell_type > internal::MatrixFreeFunctions::affine);
  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
  if (hessians_on_general_cells)
    evaluation_flag_actual |= EvaluationFlags::gradients;
 
  if (this->data->element_type ==
        internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
      evaluation_flag & EvaluationFlags::gradients &&
      (this->cell_type > internal::MatrixFreeFunctions::affine))
    evaluation_flag_actual |= EvaluationFlags::values;
 
  if constexpr (fe_degree > -1)
    internal::FEFaceEvaluationImplEvaluateInFaceSelector<
      dim,
      VectorizedArrayType>::template run<fe_degree>(n_components,
                                                    evaluation_flag_actual,
                                                    *this);
  else
    internal::FEFaceEvaluationFactory<dim, VectorizedArrayType>::
      evaluate_in_face(n_components, evaluation_flag_actual, *this);
 
  if constexpr (running_in_debug_mode())
    {
      this->values_quad_initialized =
        evaluation_flag_actual & EvaluationFlags::values;
      this->gradients_quad_initialized =
        evaluation_flag_actual & EvaluationFlags::gradients;
      this->hessians_quad_initialized =
        evaluation_flag_actual & EvaluationFlags::hessians;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
            const bool                             sum_into_values)
{
  integrate(integration_flag, this->values_dofs, sum_into_values);
 
  if constexpr (running_in_debug_mode())
    {
      this->dof_values_initialized = true;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
            VectorizedArrayType                   *values_array,
            const bool                             sum_into_values)
{
  Assert((integration_flag &
          ~(EvaluationFlags::values | EvaluationFlags::gradients |
            EvaluationFlags::hessians)) == 0,
         ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
                    "and EvaluationFlags::hessians are supported."));
 
  EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
  if (integration_flag & EvaluationFlags::hessians &&
      (this->cell_type > internal::MatrixFreeFunctions::affine))
    {
      unsigned int size = n_components * dim * n_q_points;
      if ((integration_flag & EvaluationFlags::gradients) != 0u)
        {
          for (unsigned int i = 0; i < size; ++i)
            this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
        }
      else
        {
          for (unsigned int i = 0; i < size; ++i)
            this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
          integration_flag_actual |= EvaluationFlags::gradients;
        }
    }
 
  if (this->data->element_type ==
        internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
      integration_flag & EvaluationFlags::gradients &&
      this->cell_type > internal::MatrixFreeFunctions::affine &&
      this->divergence_is_requested == false)
    {
      unsigned int size = n_components * n_q_points;
      if ((integration_flag & EvaluationFlags::values) != 0u)
        {
          for (unsigned int i = 0; i < size; ++i)
            this->values_quad[i] += this->values_from_gradients_quad[i];
        }
      else
        {
          for (unsigned int i = 0; i < size; ++i)
            this->values_quad[i] = this->values_from_gradients_quad[i];
          integration_flag_actual |= EvaluationFlags::values;
        }
    }
 
  if constexpr (fe_degree > -1)
    internal::FEFaceEvaluationImplIntegrateSelector<dim, VectorizedArrayType>::
      template run<fe_degree, n_q_points_1d>(n_components,
                                             integration_flag_actual,
                                             values_array,
                                             *this,
                                             sum_into_values);
  else
    internal::FEFaceEvaluationFactory<dim, VectorizedArrayType>::integrate(
      n_components,
      integration_flag_actual,
      values_array,
      *this,
      sum_into_values);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  integrate_in_face(const EvaluationFlags::EvaluationFlags integration_flag)
{
  Assert((integration_flag &
          ~(EvaluationFlags::values | EvaluationFlags::gradients |
            EvaluationFlags::hessians)) == 0,
         ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
                    "and EvaluationFlags::hessians are supported."));
 
  EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
  if (integration_flag & EvaluationFlags::hessians &&
      (this->cell_type > internal::MatrixFreeFunctions::affine))
    {
      unsigned int size = n_components * dim * n_q_points;
      if ((integration_flag & EvaluationFlags::gradients) != 0u)
        {
          for (unsigned int i = 0; i < size; ++i)
            this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
        }
      else
        {
          for (unsigned int i = 0; i < size; ++i)
            this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
          integration_flag_actual |= EvaluationFlags::gradients;
        }
    }
 
  if (this->data->element_type ==
        internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
      integration_flag & EvaluationFlags::gradients &&
      this->cell_type > internal::MatrixFreeFunctions::affine &&
      this->divergence_is_requested == false)
    {
      unsigned int size = n_components * n_q_points;
      if ((integration_flag & EvaluationFlags::values) != 0u)
        {
          for (unsigned int i = 0; i < size; ++i)
            this->values_quad[i] += this->values_from_gradients_quad[i];
        }
      else
        {
          for (unsigned int i = 0; i < size; ++i)
            this->values_quad[i] = this->values_from_gradients_quad[i];
          integration_flag_actual |= EvaluationFlags::values;
        }
    }
 
  if constexpr (fe_degree > -1)
    internal::FEFaceEvaluationImplIntegrateInFaceSelector<
      dim,
      VectorizedArrayType>::template run<fe_degree>(n_components,
                                                    integration_flag_actual,
                                                    *this);
  else
    internal::FEFaceEvaluationFactory<dim, VectorizedArrayType>::
      integrate_in_face(n_components, integration_flag_actual, *this);
 
  // face dofs initialized
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag,
                    const bool                             sum_into_values)
{
  collect_from_face(integration_flag, this->values_dofs, sum_into_values);
 
  if constexpr (running_in_debug_mode())
    {
      this->dof_values_initialized = true;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag,
                    VectorizedArrayType                   *values_array,
                    const bool                             sum_into_values)
{
  Assert((integration_flag &
          ~(EvaluationFlags::values | EvaluationFlags::gradients |
            EvaluationFlags::hessians)) == 0,
         ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
                    "and EvaluationFlags::hessians are supported."));
 
  EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
  if (integration_flag & EvaluationFlags::hessians &&
      (this->cell_type > internal::MatrixFreeFunctions::affine))
    integration_flag_actual |= EvaluationFlags::gradients;
 
  if (this->data->element_type ==
        internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
      integration_flag & EvaluationFlags::gradients &&
      this->cell_type > internal::MatrixFreeFunctions::affine &&
      this->divergence_is_requested == false)
    integration_flag_actual |= EvaluationFlags::values;
 
  if constexpr (fe_degree > -1)
    internal::FEFaceEvaluationImplCollectFromFaceSelector<
      dim,
      VectorizedArrayType>::template run<fe_degree>(n_components,
                                                    integration_flag_actual,
                                                    values_array,
                                                    *this,
                                                    sum_into_values);
  else
    internal::FEFaceEvaluationFactory<dim, VectorizedArrayType>::
      collect_from_face(n_components,
                        integration_flag_actual,
                        values_array,
                        *this,
                        sum_into_values);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  gather_evaluate(const VectorType                      &input_vector,
                  const EvaluationFlags::EvaluationFlags evaluation_flag)
{
  Assert((evaluation_flag &
          ~(EvaluationFlags::values | EvaluationFlags::gradients |
            EvaluationFlags::hessians)) == 0,
         ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
                    "and EvaluationFlags::hessians are supported."));
 
  const auto shared_vector_data = internal::get_shared_vector_data(
    &input_vector,
    this->dof_access_index ==
      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
    this->active_fe_index,
    this->dof_info);
 
  if (this->data->data.front().fe_degree > 0 &&
      fast_evaluation_supported(this->data->data.front().fe_degree,
                                this->data->data.front().n_q_points_1d) &&
      internal::FEFaceEvaluationImplGatherEvaluateSelector<
        dim,
        typename VectorType::value_type,
        VectorizedArrayType>::
        supports(evaluation_flag,
                 *this->data,
                 internal::get_beginning<typename VectorType::value_type>(
                   input_vector),
                 this->dof_info->index_storage_variants[this->dof_access_index]
                                                       [this->cell]))
    {
      if constexpr (fe_degree > -1)
        {
          internal::FEFaceEvaluationImplGatherEvaluateSelector<
            dim,
            typename VectorType::value_type,
            VectorizedArrayType>::template run<fe_degree,
                                               n_q_points_1d>(
            n_components,
            evaluation_flag,
            internal::get_beginning<typename VectorType::value_type>(
              input_vector),
            shared_vector_data,
            *this);
        }
      else
        {
          internal::FEFaceEvaluationGatherFactory<
            dim,
            typename VectorType::value_type,
            VectorizedArrayType>::evaluate(n_components,
                                           evaluation_flag,
                                           internal::get_beginning<
                                             typename VectorType::value_type>(
                                             input_vector),
                                           shared_vector_data,
                                           *this);
        }
    }
  else
    {
      this->read_dof_values(input_vector);
      this->evaluate(evaluation_flag);
    }
 
  if constexpr (running_in_debug_mode())
    {
      this->values_quad_initialized = evaluation_flag & EvaluationFlags::values;
      this->gradients_quad_initialized =
        evaluation_flag & EvaluationFlags::gradients;
      this->hessians_quad_initialized =
        evaluation_flag & EvaluationFlags::hessians;
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEFaceEvaluation<
  dim,
  fe_degree,
  n_q_points_1d,
  n_components_,
  Number,
  VectorizedArrayType>::integrate_scatter(const bool  integrate_values,
                                          const bool  integrate_gradients,
                                          VectorType &destination)
{
  const EvaluationFlags::EvaluationFlags flag =
    ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
    ((integrate_gradients) ? EvaluationFlags::gradients :
                             EvaluationFlags::nothing);
 
  integrate_scatter(flag, destination);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
template <typename VectorType>
inline void
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
                    VectorType                            &destination)
{
  Assert((this->dof_access_index ==
            internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
          this->is_interior_face() == false) == false,
         ExcNotImplemented());
 
  const auto shared_vector_data = internal::get_shared_vector_data(
    &destination,
    this->dof_access_index ==
      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
    this->active_fe_index,
    this->dof_info);
 
  if (this->data->data.front().fe_degree > 0 &&
      fast_evaluation_supported(this->data->data.front().fe_degree,
                                this->data->data.front().n_q_points_1d) &&
      internal::FEFaceEvaluationImplGatherEvaluateSelector<
        dim,
        typename VectorType::value_type,
        VectorizedArrayType>::
        supports(integration_flag,
                 *this->data,
                 internal::get_beginning<typename VectorType::value_type>(
                   destination),
                 this->dof_info->index_storage_variants[this->dof_access_index]
                                                       [this->cell]))
    {
      if constexpr (fe_degree > -1)
        {
          internal::FEFaceEvaluationImplIntegrateScatterSelector<
            dim,
            typename VectorType::value_type,
            VectorizedArrayType>::template run<fe_degree,
                                               n_q_points_1d>(
            n_components,
            integration_flag,
            internal::get_beginning<typename VectorType::value_type>(
              destination),
            shared_vector_data,
            *this);
        }
      else
        {
          internal::FEFaceEvaluationGatherFactory<
            dim,
            typename VectorType::value_type,
            VectorizedArrayType>::integrate(n_components,
                                            integration_flag,
                                            internal::get_beginning<
                                              typename VectorType::value_type>(
                                              destination),
                                            shared_vector_data,
                                            *this);
        }
    }
  else
    {
      integrate(integration_flag);
      this->distribute_local_to_global(destination);
    }
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::dof_indices() const
{
  return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
    0U, dofs_per_cell);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
bool
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  fast_evaluation_supported(const unsigned int given_degree,
                            const unsigned int given_n_q_points_1d)
{
  return fe_degree == -1 ?
           internal::FEEvaluationFactory<dim, VectorizedArrayType>::
             fast_evaluation_supported(given_degree, given_n_q_points_1d) :
           true;
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
bool
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::
  fast_evaluation_supported(const unsigned int given_degree,
                            const unsigned int given_n_q_points_1d)
{
  return fe_degree == -1 ?
           internal::FEFaceEvaluationFactory<dim, VectorizedArrayType>::
             fast_evaluation_supported(given_degree, given_n_q_points_1d) :
           true;
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
bool
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::at_boundary() const
{
  Assert(this->dof_access_index !=
           internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
         ExcNotImplemented());
 
  if (this->is_interior_face() == false)
    return false;
  else if (this->cell < this->matrix_free->n_inner_face_batches())
    return false;
  else if (this->cell < (this->matrix_free->n_inner_face_batches() +
                         this->matrix_free->n_boundary_face_batches()))
    return true;
  else
    return false;
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
types::boundary_id
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::boundary_id() const
{
  Assert(this->dof_access_index !=
           internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
         ExcNotImplemented());
 
  if (at_boundary())
    return this->matrix_free->get_boundary_id(this->cell);
  else
    return numbers::internal_face_boundary_id;
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
unsigned int
FEFaceEvaluation<
  dim,
  fe_degree,
  n_q_points_1d,
  n_components_,
  Number,
  VectorizedArrayType>::get_dofs_per_component_projected_to_face()
{
  return this->data->dofs_per_component_on_face;
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
unsigned int
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::get_dofs_projected_to_face()
{
  return this->data->dofs_per_component_on_face * n_components_;
}
 
 
 
/*------------------------- end FEFaceEvaluation ------------------------- */
 
 
#endif // ifndef DOXYGEN
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif