fe_evaluation.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
 
#ifndef dealii_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/smartpointer.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/vectorization.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_selector.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  = Tensor<1, n_components_, VectorizedArrayType>;
  using gradient_type =
    Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>;
  using hessian_type =
    Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>;
  static constexpr unsigned int dimension    = dim;
  static constexpr unsigned int n_components = n_components_;
 
  template <typename VectorType>
  void
  read_dof_values(const VectorType & src,
                  const unsigned int first_index = 0,
                  const std::bitset<VectorizedArrayType::size()> &mask =
                    std::bitset<VectorizedArrayType::size()>().flip());
 
  template <typename VectorType>
  void
  read_dof_values_plain(const VectorType & src,
                        const unsigned int first_index = 0,
                        const std::bitset<VectorizedArrayType::size()> &mask =
                          std::bitset<VectorizedArrayType::size()>().flip());
 
  template <typename VectorType>
  void
  distribute_local_to_global(
    VectorType &                                    dst,
    const unsigned int                              first_index = 0,
    const std::bitset<VectorizedArrayType::size()> &mask =
      std::bitset<VectorizedArrayType::size()>().flip()) const;
 
  template <typename VectorType>
  void
  set_dof_values(VectorType &       dst,
                 const unsigned int first_index = 0,
                 const std::bitset<VectorizedArrayType::size()> &mask =
                   std::bitset<VectorizedArrayType::size()>().flip()) const;
 
  template <typename VectorType>
  void
  set_dof_values_plain(
    VectorType &                                    dst,
    const unsigned int                              first_index = 0,
    const std::bitset<VectorizedArrayType::size()> &mask =
      std::bitset<VectorizedArrayType::size()>().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);
 
  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);
 
  void
  submit_normal_derivative(const value_type   grad_in,
                           const unsigned int q_point);
 
  void
  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
 
  Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
  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;
 
#ifdef DOXYGEN
  // doxygen does not anyhow mention functions coming from partial template
  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
  // For now, hack in those functions manually only to fix documentation:
 
  VectorizedArrayType
  get_divergence(const unsigned int q_point) const;
 
  SymmetricTensor<2, dim, VectorizedArrayType>
  get_symmetric_gradient(const unsigned int q_point) const;
 
  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
  get_curl(const unsigned int q_point) const;
 
  void
  submit_divergence(const VectorizedArrayType div_in,
                    const unsigned int        q_point);
 
  void
  submit_symmetric_gradient(
    const SymmetricTensor<2, dim, VectorizedArrayType> grad_in,
    const unsigned int                                 q_point);
 
  void
  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
              const unsigned int                                       q_point);
 
#endif
 
  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_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no,
    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<VectorizedArrayType::size()> &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<VectorizedArrayType::size()> &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 n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType = VectorizedArray<Number>>
class FEEvaluationAccess : public FEEvaluationBase<dim,
                                                   n_components_,
                                                   Number,
                                                   is_face,
                                                   VectorizedArrayType>
{
  static_assert(
    std::is_same<Number, typename VectorizedArrayType::value_type>::value,
    "Type of Number and of VectorizedArrayType do not match.");
 
public:
  using number_type = Number;
  using value_type  = Tensor<1, n_components_, VectorizedArrayType>;
  using gradient_type =
    Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>;
  static constexpr unsigned int dimension    = dim;
  static constexpr unsigned int n_components = n_components_;
  using BaseClass =
    FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>;
 
protected:
  FEEvaluationAccess(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                  dof_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no,
    const unsigned int fe_degree,
    const unsigned int n_q_points,
    const bool         is_interior_face  = true,
    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);
 
  FEEvaluationAccess(
    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);
 
  FEEvaluationAccess(const FEEvaluationAccess &other);
 
  FEEvaluationAccess &
  operator=(const FEEvaluationAccess &other);
};
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
  : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
{
  static_assert(
    std::is_same<Number, typename VectorizedArrayType::value_type>::value,
    "Type of Number and of VectorizedArrayType do not match.");
 
public:
  using number_type                       = Number;
  using value_type                        = VectorizedArrayType;
  using gradient_type                     = Tensor<1, dim, VectorizedArrayType>;
  using hessian_type                      = Tensor<2, dim, VectorizedArrayType>;
  static constexpr unsigned int dimension = dim;
  using BaseClass =
    FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>;
 
  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);
 
  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);
 
  void
  submit_normal_derivative(const value_type   grad_in,
                           const unsigned int q_point);
 
  hessian_type
  get_hessian(unsigned int q_point) const;
 
  gradient_type
  get_hessian_diagonal(const unsigned int q_point) const;
 
  void
  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
 
  value_type
  get_laplacian(const unsigned int q_point) const;
 
  value_type
  integrate_value() const;
 
protected:
  FEEvaluationAccess(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                  dof_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no,
    const unsigned int fe_degree,
    const unsigned int n_q_points,
    const bool         is_interior_face  = true,
    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);
 
  FEEvaluationAccess(
    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);
 
  FEEvaluationAccess(const FEEvaluationAccess &other);
 
  FEEvaluationAccess &
  operator=(const FEEvaluationAccess &other);
};
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
  : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
{
  static_assert(
    std::is_same<Number, typename VectorizedArrayType::value_type>::value,
    "Type of Number and of VectorizedArrayType do not match.");
 
public:
  using number_type                       = Number;
  using value_type                        = Tensor<1, dim, VectorizedArrayType>;
  using gradient_type                     = Tensor<2, dim, VectorizedArrayType>;
  static constexpr unsigned int dimension = dim;
  static constexpr unsigned int n_components = dim;
  using BaseClass =
    FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>;
 
  value_type
  get_value(const unsigned int q_point) const;
 
  gradient_type
  get_gradient(const unsigned int q_point) const;
 
  VectorizedArrayType
  get_divergence(const unsigned int q_point) const;
 
  SymmetricTensor<2, dim, VectorizedArrayType>
  get_symmetric_gradient(const unsigned int q_point) const;
 
  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
  get_curl(const unsigned int q_point) const;
 
  Tensor<3, dim, VectorizedArrayType>
  get_hessian(const unsigned int q_point) const;
 
  gradient_type
  get_hessian_diagonal(const unsigned int q_point) const;
 
  void
  submit_value(const Tensor<1, dim, VectorizedArrayType> val_in,
               const unsigned int                        q_point);
 
  void
  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
 
  void
  submit_gradient(
    const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
    const unsigned int                                        q_point);
 
  void
  submit_divergence(const VectorizedArrayType div_in,
                    const unsigned int        q_point);
 
  void
  submit_symmetric_gradient(
    const SymmetricTensor<2, dim, VectorizedArrayType> grad_in,
    const unsigned int                                 q_point);
 
  void
  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
              const unsigned int                                       q_point);
 
protected:
  FEEvaluationAccess(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                  dof_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no,
    const unsigned int dofs_per_cell,
    const unsigned int n_q_points,
    const bool         is_interior_face  = true,
    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);
 
  FEEvaluationAccess(
    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);
 
  FEEvaluationAccess(const FEEvaluationAccess &other);
 
  FEEvaluationAccess &
  operator=(const FEEvaluationAccess &other);
};
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
  : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
{
  static_assert(
    std::is_same<Number, typename VectorizedArrayType::value_type>::value,
    "Type of Number and of VectorizedArrayType do not match.");
 
public:
  using number_type                       = Number;
  using value_type                        = VectorizedArrayType;
  using gradient_type                     = Tensor<1, 1, VectorizedArrayType>;
  using hessian_type                      = Tensor<2, 1, VectorizedArrayType>;
  static constexpr unsigned int dimension = 1;
  using BaseClass =
    FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>;
 
  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);
 
  void
  submit_value(const gradient_type val_in, const unsigned int q_point);
 
  gradient_type
  get_gradient(const unsigned int q_point) const;
 
  value_type
  get_divergence(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);
 
  void
  submit_gradient(const value_type grad_in, const unsigned int q_point);
 
  void
  submit_gradient(const Tensor<2, 1, VectorizedArrayType> grad_in,
                  const unsigned int                      q_point);
 
  void
  submit_normal_derivative(const value_type   grad_in,
                           const unsigned int q_point);
 
  void
  submit_normal_derivative(const gradient_type grad_in,
                           const unsigned int  q_point);
 
  hessian_type
  get_hessian(unsigned int q_point) const;
 
  gradient_type
  get_hessian_diagonal(const unsigned int q_point) const;
 
  void
  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
 
  value_type
  get_laplacian(const unsigned int q_point) const;
 
  value_type
  integrate_value() const;
 
protected:
  FEEvaluationAccess(
    const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                dof_no,
    const unsigned int                                first_selected_component,
    const unsigned int                                quad_no,
    const unsigned int                                fe_degree,
    const unsigned int                                n_q_points,
    const bool                                        is_interior_face = true,
    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);
 
  FEEvaluationAccess(
    const Mapping<1> &      mapping,
    const FiniteElement<1> &fe,
    const Quadrature<1> &   quadrature,
    const UpdateFlags       update_flags,
    const unsigned int      first_selected_component,
    const FEEvaluationData<1, VectorizedArrayType, is_face> *other);
 
  FEEvaluationAccess(const FEEvaluationAccess &other);
 
  FEEvaluationAccess &
  operator=(const FEEvaluationAccess &other);
};
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
class FEEvaluation : public FEEvaluationAccess<dim,
                                               n_components_,
                                               Number,
                                               false,
                                               VectorizedArrayType>
{
  static_assert(
    std::is_same<Number, typename VectorizedArrayType::value_type>::value,
    "Type of Number and of VectorizedArrayType do not match.");
 
public:
  using BaseClass =
    FEEvaluationAccess<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 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_no  = 0,
    const unsigned int                                  quad_no = 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_no  = 0,
               const unsigned int                                  quad_no = 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, VectorizedArrayType::size()> &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);
 
  DEAL_II_DEPRECATED void
  evaluate(const bool evaluate_values,
           const bool evaluate_gradients,
           const bool evaluate_hessians = false);
 
  void
  evaluate(const VectorizedArrayType *            values_array,
           const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  DEAL_II_DEPRECATED void
  evaluate(const VectorizedArrayType *values_array,
           const bool                 evaluate_values,
           const bool                 evaluate_gradients,
           const bool                 evaluate_hessians = false);
 
  template <typename VectorType>
  void
  gather_evaluate(const VectorType &                     input_vector,
                  const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  template <typename VectorType>
  DEAL_II_DEPRECATED void
  gather_evaluate(const VectorType &input_vector,
                  const bool        evaluate_values,
                  const bool        evaluate_gradients,
                  const bool        evaluate_hessians = false);
 
  void
  integrate(const EvaluationFlags::EvaluationFlags integration_flag);
 
  DEAL_II_DEPRECATED void
  integrate(const bool integrate_values, const bool integrate_gradients);
 
  void
  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
            VectorizedArrayType *                  values_array,
            const bool                             sum_into_values = false);
 
  DEAL_II_DEPRECATED void
  integrate(const bool           integrate_values,
            const bool           integrate_gradients,
            VectorizedArrayType *values_array);
 
  template <typename VectorType>
  void
  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
                    VectorType &                           output_vector);
 
  template <typename VectorType>
  DEAL_II_DEPRECATED 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;
 
  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 FEEvaluationAccess<dim,
                                                   n_components_,
                                                   Number,
                                                   true,
                                                   VectorizedArrayType>
{
  static_assert(
    std::is_same<Number, typename VectorizedArrayType::value_type>::value,
    "Type of Number and of VectorizedArrayType do not match.");
 
public:
  using BaseClass =
    FEEvaluationAccess<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 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_no           = 0,
    const unsigned int                                  quad_no          = 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_no           = 0,
    const unsigned int                                  quad_no          = 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);
 
  DEAL_II_DEPRECATED void
  evaluate(const bool evaluate_values, const bool evaluate_gradients);
 
  void
  evaluate(const VectorizedArrayType *            values_array,
           const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  DEAL_II_DEPRECATED void
  evaluate(const VectorizedArrayType *values_array,
           const bool                 evaluate_values,
           const bool                 evaluate_gradients);
 
  template <typename VectorType>
  void
  gather_evaluate(const VectorType &                     input_vector,
                  const EvaluationFlags::EvaluationFlags evaluation_flag);
 
  template <typename VectorType>
  DEAL_II_DEPRECATED void
  gather_evaluate(const VectorType &input_vector,
                  const bool        evaluate_values,
                  const bool        evaluate_gradients);
 
  void
  integrate(const EvaluationFlags::EvaluationFlags integration_flag);
 
  DEAL_II_DEPRECATED void
  integrate(const bool integrate_values, const bool integrate_gradients);
 
  void
  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
            VectorizedArrayType *                  values_array);
 
  DEAL_II_DEPRECATED void
  integrate(const bool           integrate_values,
            const bool           integrate_gradients,
            VectorizedArrayType *values_array);
 
  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;
 
  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_no,
      const unsigned int first_selected_component,
      const unsigned int quad_no,
      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_no);
    init_data.mapping_data =
      &internal::MatrixFreeFunctions::
        MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
          matrix_free.get_mapping_info(), quad_no);
 
    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() -
                                    1)) :
        init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
 
    init_data.shape_info = &matrix_free.get_shape_info(
      dof_no,
      quad_no,
      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_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no,
    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_no,
                                                     first_selected_component,
                                                     quad_no,
                                                     fe_degree,
                                                     n_q_points,
                                                     active_fe_index,
                                                     active_quad_index,
                                                     face_type),
      is_interior_face,
      quad_no,
      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<VectorizedArrayType>(
      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<VectorizedArrayType>(
          *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());
      this->mapping_data = &this->mapped_geometry->get_data_storage();
      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<VectorizedArrayType>(
          *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());
      this->cell         = 0;
      this->mapping_data = &this->mapped_geometry->get_data_storage();
      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<VectorType>::value>::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];
    }
  };
} // 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<VectorizedArrayType::size()> &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);
    }
 
  // 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 (is_face && !this->interior_face &&
          (this->dof_access_index ==
           internal::MatrixFreeFunctions::DoFInfo::dof_access_cell))
        {
          const std::array<unsigned int, VectorizedArrayType::size()> &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] / VectorizedArrayType::size()] <
                  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
  constexpr unsigned int n_lanes = VectorizedArrayType::size();
 
  std::array<unsigned int, VectorizedArrayType::size()> 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;
 
  bool has_hn_constraints = false;
 
  if (is_face == false)
    {
      if (!dof_info.hanging_node_constraint_masks.empty() &&
          !dof_info.hanging_node_constraint_masks_comp.empty() &&
          dof_info
            .hanging_node_constraint_masks_comp[this->active_fe_index]
                                               [this->first_selected_component])
        for (unsigned int v = 0; v < n_lanes; ++v)
          if (cells[v] != numbers::invalid_unsigned_int &&
              dof_info.hanging_node_constraint_masks[cells[v]] !=
                internal::MatrixFreeFunctions::
                  unconstrained_compressed_constraint_kind)
            has_hn_constraints = true;
    }
 
  std::integral_constant<bool,
                         internal::is_vectorizable<VectorType, Number>::value>
    vector_selector;
 
  const bool use_vectorized_path = !(masking_is_active || has_hn_constraints);
 
  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
  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;
 
  if (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;
      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,
                                         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, dof_indices += n_lanes)
            operation.process_dof_gather(
              dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
      return;
    }
 
  // 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;
 
          if (dof_info.hanging_node_constraint_masks.size() > 0 &&
              dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
              dof_info.hanging_node_constraint_masks[cells[v]] !=
                internal::MatrixFreeFunctions::
                  unconstrained_compressed_constraint_kind &&
              dof_info.hanging_node_constraint_masks_comp
                [this->active_fe_index][this->first_selected_component])
            has_hn_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;
        }
    }
 
  if (std::count_if(cells.begin(), cells.end(), [](const auto i) {
        return i != numbers::invalid_unsigned_int;
      }) < n_lanes)
    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)
                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)
                  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;
      const unsigned int n_components_read =
        this->n_fe_components > 1 ? n_components : 1;
      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[cell_dof_index].second !=
             dof_info.row_starts[cell_dof_index + n_components_read].second ||
           ((dof_info.hanging_node_constraint_masks.size() > 0 &&
             dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
             dof_info.hanging_node_constraint_masks[cell_index] !=
               internal::MatrixFreeFunctions::
                 unconstrained_compressed_constraint_kind) &&
            dof_info.hanging_node_constraint_masks_comp
              [this->active_fe_index][this->first_selected_component])))
        {
          Assert(dof_info.row_starts_plain_indices[cell_index] !=
                   numbers::invalid_unsigned_int,
                 ExcNotInitialized());
          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)
                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])
                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)
            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)
                    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])
                    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)
                {
                  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<VectorizedArrayType::size()> &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::integral_constant<bool,
                         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_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;
  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());
 
  // 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_lanes == VectorizedArrayType::size()) &&
      !(is_face &&
        this->dof_access_index ==
          internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
        this->is_interior_face() == false) &&
      !(!is_face && !this->is_interior_face()))
    {
      const unsigned int dof_index =
        dof_indices_cont[this->cell * VectorizedArrayType::size()] +
        this->dof_info
            ->component_dof_indices_offset[this->active_fe_index]
                                          [this->first_selected_component] *
          VectorizedArrayType::size();
      if (n_components == 1 || this->n_fe_components == 1)
        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, VectorizedArrayType::size()> &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 is_ecl =
    (this->dof_access_index ==
       internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
     this->is_interior_face() == false) ||
    (!is_face && !this->is_interior_face());
 
  if (vectors_sm[0] != nullptr)
    {
      const auto compute_vector_ptrs = [&](const unsigned int comp) {
        std::array<typename VectorType::value_type *,
                   VectorizedArrayType::size()>
          vector_ptrs = {};
 
        const auto upper_bound =
          std::min<unsigned int>(n_filled_lanes, VectorizedArrayType::size());
        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 < VectorizedArrayType::size();
             ++v)
          vector_ptrs[v] = nullptr;
 
        return vector_ptrs;
      };
 
      if (n_filled_lanes == VectorizedArrayType::size() &&
          n_lanes == VectorizedArrayType::size() && !is_ecl)
        {
          if (n_components == 1 || this->n_fe_components == 1)
            {
              for (unsigned int comp = 0; comp < n_components; ++comp)
                {
                  auto vector_ptrs = compute_vector_ptrs(comp);
                  operation.process_dofs_vectorized_transpose(
                    dofs_per_component,
                    vector_ptrs,
                    values_dofs[comp],
                    vector_selector);
                }
            }
          else
            {
              auto vector_ptrs = compute_vector_ptrs(0);
              operation.process_dofs_vectorized_transpose(dofs_per_component *
                                                            n_components,
                                                          vector_ptrs,
                                                          &values_dofs[0][0],
                                                          vector_selector);
            }
        }
      else
        for (unsigned int comp = 0; comp < n_components; ++comp)
          {
            auto vector_ptrs = compute_vector_ptrs(
              (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
 
            for (unsigned int i = 0; i < dofs_per_component; ++i)
              operation.process_empty(values_dofs[comp][i]);
 
            if (n_components == 1 || this->n_fe_components == 1)
              {
                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],
                                            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(
                        vector_ptrs[v][i + comp * dofs_per_component],
                        values_dofs[comp][i][v]);
              }
          }
      return;
    }
 
  unsigned int dof_indices[VectorizedArrayType::size()];
 
  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]];
      else
        dof_indices[v] = numbers::invalid_unsigned_int;
    }
 
  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
    dof_indices[v] = numbers::invalid_unsigned_int;
 
  // In the case with contiguous cell indices, we know that there are no
  // constraints and that the indices within each element are contiguous
  if (n_filled_lanes == VectorizedArrayType::size() &&
      n_lanes == VectorizedArrayType::size() && !is_ecl)
    {
      if (dof_info.index_storage_variants[ind][this->cell] ==
          internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
            contiguous)
        {
          if (n_components == 1 || this->n_fe_components == 1)
            for (unsigned int comp = 0; comp < n_components; ++comp)
              operation.process_dofs_vectorized_transpose(dofs_per_component,
                                                          dof_indices,
                                                          *src[comp],
                                                          values_dofs[comp],
                                                          vector_selector);
          else
            operation.process_dofs_vectorized_transpose(dofs_per_component *
                                                          n_components,
                                                        dof_indices,
                                                        *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)
        {
          if (n_components == 1 || this->n_fe_components == 1)
            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,
                                               *src[comp],
                                               i * VectorizedArrayType::size(),
                                               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,
                                               *src[0],
                                               (comp * dofs_per_component + i) *
                                                 VectorizedArrayType::size(),
                                               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());
          const unsigned int *offsets =
            &dof_info.dof_indices_interleave_strides
               [ind][VectorizedArrayType::size() * this->cell];
          if (n_components == 1 || this->n_fe_components == 1)
            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,
                                               *src[comp],
                                               0,
                                               values_dofs[comp][i],
                                               vector_selector);
                DEAL_II_OPENMP_SIMD_PRAGMA
                for (unsigned int v = 0; v < VectorizedArrayType::size(); ++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,
                                               *src[0],
                                               0,
                                               values_dofs[comp][i],
                                               vector_selector);
                  DEAL_II_OPENMP_SIMD_PRAGMA
                  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++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 (dof_info.index_storage_variants[ind][this->cell] ==
            internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
              contiguous)
          {
            if (n_components == 1 || this->n_fe_components == 1)
              {
                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 (n_components == 1 || this->n_fe_components == 1)
              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<
          typename std::remove_const<VectorType>::type,
          IsBlockVector<typename std::remove_const<VectorType>::type>::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.size() == 0 ||
      this->dof_info->hanging_node_constraint_masks_comp.size() == 0 ||
      this->dof_info->hanging_node_constraint_masks_comp
          [this->active_fe_index][this->first_selected_component] == false)
    return; // nothing to do with faces
 
  constexpr unsigned int n_lanes = VectorizedArrayType::size();
  std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
    constraint_mask;
 
  bool hn_available = false;
 
  const std::array<unsigned int, VectorizedArrayType::size()> &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<VectorizedArrayType::size()> &mask)
{
  const auto src_data = internal::get_vector_data<n_components_>(
    src,
    first_index,
    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);
 
#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif
}
 
 
 
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<VectorizedArrayType::size()> &mask)
{
  const auto src_data = internal::get_vector_data<n_components_>(
    src,
    first_index,
    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);
 
#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif
}
 
 
 
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<VectorizedArrayType::size()> &mask) const
{
#  ifdef DEBUG
  Assert(this->dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
 
  apply_hanging_node_constraints(true);
 
  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::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<VectorizedArrayType::size()> &mask) const
{
#  ifdef DEBUG
  Assert(this->dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
 
  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);
}
 
 
 
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<VectorizedArrayType::size()> &mask) const
{
#  ifdef DEBUG
  Assert(this->dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
 
  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 Tensor<1, n_components_, VectorizedArrayType>
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);
  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 Tensor<1, n_components_, VectorizedArrayType>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_value(const unsigned int q_point) const
{
#  ifdef DEBUG
  Assert(this->values_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
 
  AssertIndexRange(q_point, this->n_quadrature_points);
  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
  Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>
  FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
    get_gradient(const unsigned int q_point) const
{
#  ifdef DEBUG
  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
 
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  const std::size_t nqp = this->n_quadrature_points;
  Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_out;
 
  // Cartesian cell
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      for (unsigned int d = 0; d < dim; ++d)
        for (unsigned int comp = 0; comp < n_components; ++comp)
          grad_out[comp][d] =
            this->gradients_quad[(comp * dim + d) * nqp + q_point] *
            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 * dim) * nqp + q_point];
            for (unsigned int e = 1; e < dim; ++e)
              grad_out[comp][d] +=
                jac[d][e] *
                this->gradients_quad[(comp * dim + e) * nqp + q_point];
          }
    }
  return grad_out;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArrayType>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_normal_derivative(const unsigned int q_point) const
{
  AssertIndexRange(q_point, this->n_quadrature_points);
#  ifdef DEBUG
  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
 
  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 * dim + dim - 1) * nqp + q_point] *
        (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 * dim * nqp + q_point] *
                           this->normal_x_jacobian[index][0];
          for (unsigned int d = 1; d < dim; ++d)
            grad_out[comp] +=
              this->gradients_quad[(comp * dim + d) * nqp + q_point] *
              this->normal_x_jacobian[index][d];
        }
    }
  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 Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_hessian(const unsigned int q_point) const
{
#  ifdef DEBUG
  Assert(this->hessians_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  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:
                Assert(false, ExcNotImplemented());
            }
          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 * dim + e) * nqp + q_point];
 
          // 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 * dim + f) * nqp + q_point];
 
          // 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];
        }
    }
  return hessian_out;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_hessian_diagonal(const unsigned int q_point) const
{
  Assert(!is_face, ExcNotImplemented());
#  ifdef DEBUG
  Assert(this->hessians_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  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 * dim + e) * nqp + q_point];
        }
    }
  return hessian_out;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline Tensor<1, n_components_, VectorizedArrayType>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_laplacian(const unsigned int q_point) const
{
  Assert(is_face == false, ExcNotImplemented());
#  ifdef DEBUG
  Assert(this->hessians_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
 
  Tensor<1, n_components_, VectorizedArrayType> laplacian_out;
  const auto hess_diag = get_hessian_diagonal(q_point);
  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 DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  submit_dof_value(const Tensor<1, n_components_, VectorizedArrayType> val_in,
                   const unsigned int                                  dof)
{
#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif
  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)
    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 Tensor<1, n_components_, VectorizedArrayType> val_in,
               const unsigned int                                  q_point)
{
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_values"));
#  ifdef DEBUG
  this->values_quad_submitted = true;
#  endif
 
  const std::size_t nqp = this->n_quadrature_points;
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
    }
  else
    {
      const VectorizedArrayType JxW = this->J_value[q_point];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        this->values_quad[comp * nqp + q_point] = val_in[comp] * 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_gradient(
    const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
    const unsigned int                                                  q_point)
{
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  const std::size_t nqp = this->n_quadrature_points;
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
      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 = jac[d] * JxW;
          for (unsigned int comp = 0; comp < n_components; ++comp)
            this->gradients_quad[(comp * dim + d) * nqp + q_point] =
              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];
      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]);
            this->gradients_quad[(comp * dim + d) * nqp + q_point] =
              new_val * 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_derivative(
    const Tensor<1, n_components_, VectorizedArrayType> grad_in,
    const unsigned int                                  q_point)
{
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->normal_x_jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  const std::size_t nqp = this->n_quadrature_points;
  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)
            this->gradients_quad[(comp * dim + d) * nqp + q_point] =
              VectorizedArrayType();
          this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
            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)
            this->gradients_quad[(comp * dim + d) * nqp + q_point] =
              (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 Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
                       hessian_in,
    const unsigned int q_point)
{
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_hessians"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_hessians"));
#  ifdef DEBUG
  this->hessians_quad_submitted = true;
#  endif
 
  // 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)
            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)
              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)
        {
          // 1. tmp = hessian_in(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_in[comp][i][0] * jac[0][j];
                for (unsigned int k = 1; k < dim; ++k)
                  tmp[i][j] += hessian_in[comp][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)
        {
          // 1. tmp = hessian_in(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_in[comp][i][0] * jac[0][j];
                for (unsigned int k = 1; k < dim; ++k)
                  tmp[i][j] += hessian_in[comp][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_in
          for (unsigned int d = 0; d < dim; ++d)
            {
              VectorizedArrayType sum = 0;
              for (unsigned int e = 0; e < dim; ++e)
                sum += hessian_in[comp][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_in[comp][e][f] + hessian_in[comp][f][e]) *
                         jac_grad[count][d];
              this->gradients_from_hessians_quad[(comp * dim + d) * nqp +
                                                 q_point] = sum * JxW;
            }
        }
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline Tensor<1, n_components_, VectorizedArrayType>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  integrate_value() const
{
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
  Assert(this->values_quad_submitted == true,
         internal::ExcAccessToUninitializedField());
#  endif
 
  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];
  return (return_value);
}
 
 
 
/*----------------------- FEEvaluationAccess --------------------------------*/
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline FEEvaluationAccess<dim,
                          n_components_,
                          Number,
                          is_face,
                          VectorizedArrayType>::
  FEEvaluationAccess(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                  dof_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no,
    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<dim, n_components_, Number, is_face, VectorizedArrayType>(
      matrix_free,
      dof_no,
      first_selected_component,
      quad_no,
      fe_degree,
      n_q_points,
      is_interior_face,
      active_fe_index,
      active_quad_index,
      face_type)
{}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline FEEvaluationAccess<dim,
                          n_components_,
                          Number,
                          is_face,
                          VectorizedArrayType>::
  FEEvaluationAccess(
    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<dim, n_components_, Number, is_face, VectorizedArrayType>(
      mapping,
      fe,
      quadrature,
      update_flags,
      first_selected_component,
      other)
{}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline FEEvaluationAccess<dim,
                          n_components_,
                          Number,
                          is_face,
                          VectorizedArrayType>::
  FEEvaluationAccess(const FEEvaluationAccess<dim,
                                              n_components_,
                                              Number,
                                              is_face,
                                              VectorizedArrayType> &other)
  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
      other)
{}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline FEEvaluationAccess<dim,
                          n_components_,
                          Number,
                          is_face,
                          VectorizedArrayType> &
FEEvaluationAccess<dim, n_components_, Number, is_face, VectorizedArrayType>::
operator=(const FEEvaluationAccess<dim,
                                   n_components_,
                                   Number,
                                   is_face,
                                   VectorizedArrayType> &other)
{
  this->FEEvaluationBase<dim,
                         n_components_,
                         Number,
                         is_face,
                         VectorizedArrayType>::operator=(other);
  return *this;
}
 
 
 
/*-------------------- FEEvaluationAccess scalar ----------------------------*/
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  FEEvaluationAccess(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                  dof_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no,
    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<dim, 1, Number, is_face, VectorizedArrayType>(
      matrix_free,
      dof_no,
      first_selected_component,
      quad_no,
      fe_degree,
      n_q_points,
      is_interior_face,
      active_fe_index,
      active_quad_index,
      face_type)
{}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  FEEvaluationAccess(
    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<dim, 1, Number, is_face, VectorizedArrayType>(
      mapping,
      fe,
      quadrature,
      update_flags,
      first_selected_component,
      other)
{}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  FEEvaluationAccess(
    const FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
      &other)
  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(other)
{}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType> &
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
  const FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType> &other)
{
  this
    ->FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
      other);
  return *this;
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::get_dof_value(
  const unsigned int dof) const
{
  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
  return this->values_dofs[dof];
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::get_value(
  const unsigned int q_point) const
{
#  ifdef DEBUG
  Assert(this->values_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  return this->values_quad[q_point];
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  get_normal_derivative(const unsigned int q_point) const
{
  return BaseClass::get_normal_derivative(q_point)[0];
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::get_gradient(
  const unsigned int q_point) const
{
  // could use the base class gradient, but that involves too many expensive
  // initialization operations on tensors
 
#  ifdef DEBUG
  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
 
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
 
  Tensor<1, dim, VectorizedArrayType> grad_out;
 
  const std::size_t nqp = this->n_quadrature_points;
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      for (unsigned int d = 0; d < dim; ++d)
        grad_out[d] =
          this->gradients_quad[d * nqp + q_point] * 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 d = 0; d < dim; ++d)
        {
          grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
          for (unsigned int e = 1; e < dim; ++e)
            grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
        }
    }
  return grad_out;
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline Tensor<2, dim, VectorizedArrayType>
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::get_hessian(
  const unsigned int q_point) const
{
  return BaseClass::get_hessian(q_point)[0];
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline Tensor<1, dim, VectorizedArrayType>
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  get_hessian_diagonal(const unsigned int q_point) const
{
  return BaseClass::get_hessian_diagonal(q_point)[0];
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline VectorizedArrayType
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::get_laplacian(
  const unsigned int q_point) const
{
  return BaseClass::get_laplacian(q_point)[0];
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline void DEAL_II_ALWAYS_INLINE
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
{
#  ifdef DEBUG
  this->dof_values_initialized = true;
  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
#  endif
  this->values_dofs[dof] = val_in;
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline void DEAL_II_ALWAYS_INLINE
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::submit_value(
  const VectorizedArrayType val_in,
  const unsigned int        q_point)
{
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_value"));
#  ifdef DEBUG
  this->values_quad_submitted = true;
#  endif
 
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
      this->values_quad[q_point] = val_in * JxW;
    }
  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
    {
      this->values_quad[q_point] = val_in * this->J_value[q_point];
    }
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::submit_value(
  const Tensor<1, 1, VectorizedArrayType> val_in,
  const unsigned int                      q_point)
{
  submit_value(val_in[0], q_point);
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  submit_normal_derivative(const VectorizedArrayType grad_in,
                           const unsigned int        q_point)
{
  Tensor<1, 1, VectorizedArrayType> grad;
  grad[0] = grad_in;
  BaseClass::submit_normal_derivative(grad, q_point);
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  submit_gradient(const Tensor<1, dim, VectorizedArrayType> grad_in,
                  const unsigned int                        q_point)
{
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  const std::size_t nqp = this->n_quadrature_points;
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
 
      // Make sure the compiler does not think 'jacobian' is aliased with
      // 'gradients_quad'
      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)
        this->gradients_quad[d * nqp + q_point] = grad_in[d] * jac[d] * JxW;
    }
  // general/affine cell type
  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];
      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];
          this->gradients_quad[d * nqp + q_point] = new_val * JxW;
        }
    }
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  submit_hessian(const Tensor<2, dim, VectorizedArrayType> hessian_in,
                 const unsigned int                        q_point)
{
  Tensor<1, 1, Tensor<2, dim, VectorizedArrayType>> hessian;
  hessian[0] = hessian_in;
  BaseClass::submit_hessian(hessian, q_point);
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline VectorizedArrayType
FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
  integrate_value() const
{
  return BaseClass::integrate_value()[0];
}
 
 
 
/*----------------- FEEvaluationAccess vector-valued ------------------------*/
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  FEEvaluationAccess(
    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                  dof_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no,
    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<dim, dim, Number, is_face, VectorizedArrayType>(
      matrix_free,
      dof_no,
      first_selected_component,
      quad_no,
      fe_degree,
      n_q_points,
      is_interior_face,
      active_fe_index,
      active_quad_index,
      face_type)
{}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  FEEvaluationAccess(
    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<dim, dim, Number, is_face, VectorizedArrayType>(
      mapping,
      fe,
      quadrature,
      update_flags,
      first_selected_component,
      other)
{}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  FEEvaluationAccess(
    const FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
      &other)
  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(other)
{}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType> &
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::operator=(
  const FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
    &other)
{
  this->FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>::
  operator=(other);
  return *this;
}
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_value(
  const unsigned int q_point) const
{
  if (this->data->element_type ==
      internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
    {
      // Piola transform is required
#  ifdef DEBUG
      Assert(this->values_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
#  endif
 
      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, dim, 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] * inv_det;
              for (unsigned int e = 1; e < dim; ++e)
                value_out[comp] +=
                  this->values_quad[e * nqp + q_point] * jac[comp][e] * inv_det;
            }
        }
      return value_out;
    }
  else
    {
      // No Piola needed
      return BaseClass::get_value(q_point);
    }
}
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  get_gradient(const unsigned int q_point) const
{
  if (this->data->element_type ==
      internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
    {
      // Piola transform is required
#  ifdef DEBUG
      Assert(this->gradients_quad_initialized == true,
             internal::ExcAccessToUninitializedField());
#  endif
 
      AssertIndexRange(q_point, this->n_quadrature_points);
      Assert(this->jacobian != nullptr,
             internal::ExcMatrixFreeAccessToUninitializedMappingField(
               "update_gradients"));
      const std::size_t nqp = this->n_quadrature_points;
      Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_out;
 
      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] =
                this->gradients_quad[(comp * dim + d) * nqp + q_point] *
                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;
          // J * grad_quad * J^-1 * det(J^-1)
          for (unsigned int comp = 0; comp < n_components; ++comp)
            for (unsigned int d = 0; d < dim; ++d)
              {
                tmp = 0;
                for (unsigned int f = 0; f < dim; ++f)
                  for (unsigned int e = 0; e < dim; ++e)
                    tmp += jac[comp][f] * inv_t_jac[d][e] * inv_det *
                           this->gradients_quad[(f * dim + e) * nqp + q_point];
 
                grad_out[comp][d] = tmp;
              }
        }
      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];
          const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
 
          // 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;
          // J * grad_quad * J^-1 * det(J^-1)
          for (unsigned int comp = 0; comp < n_components; ++comp)
            for (unsigned int d = 0; d < dim; ++d)
              {
                tmp = 0;
                for (unsigned int f = 0; f < dim; ++f)
                  for (unsigned int e = 0; e < dim; ++e)
                    tmp += t_jac[f][comp] * inv_t_jac[d][e] *
                           this->gradients_quad[(f * dim + e) * nqp + q_point];
 
                grad_out[comp][d] = tmp * inv_det;
              }
 
          // Contribution from values
          {
            // Diagonal part of jac_grad
 
            // Add jac_grad * J^{-1} * values * det(J^{-1})
            // -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1}))
            for (unsigned int i = 0; i < dim; ++i)
              for (unsigned int j = 0; j < dim; ++j)
                {
                  tmp = jac_grad[0][i] * inv_t_jac[j][0] *
                        this->values_quad[q_point];
                  for (unsigned int f = 1; f < dim; ++f)
                    tmp += jac_grad[f][i] * inv_t_jac[j][f] *
                           this->values_quad[f * nqp + q_point];
 
                  grad_out[i][j] += tmp * inv_det;
                }
 
            for (unsigned int i = 0; i < dim; ++i)
              for (unsigned int j = 0; j < dim; ++j)
                {
                  tmp = 0;
                  for (unsigned int f = 0; f < dim; ++f)
                    for (unsigned int n = 0; n < dim; ++n)
                      for (unsigned int m = 0; m < dim; ++m)
                        tmp += inv_t_jac[m][f] * jac_grad[f][m] *
                               inv_t_jac[j][f] * t_jac[n][i] *
                               this->values_quad[n * nqp + q_point];
                  grad_out[i][j] -= tmp * inv_det;
                }
          }
 
          {
            // Off-diagonal part of jac_grad
 
            // Add jac_grad * J^{-1} * values * det(J^{-1})
            // -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1}))
            for (unsigned int i = 0; i < dim; ++i)
              for (unsigned int j = 0; j < dim; ++j)
                {
                  tmp = 0;
                  for (unsigned int r = 0, f = dim; r < dim; ++r)
                    for (unsigned int k = r + 1; k < dim; ++k, ++f)
                      {
                        tmp += jac_grad[f][i] *
                               (inv_t_jac[j][k] *
                                  this->values_quad[r * nqp + q_point] +
                                inv_t_jac[j][r] *
                                  this->values_quad[k * nqp + q_point]);
                        for (unsigned int n = 0; n < dim; ++n)
                          for (unsigned int m = 0; m < dim; ++m)
                            tmp -= jac_grad[f][m] * t_jac[n][i] *
                                   this->values_quad[n * nqp + q_point] *
                                   (inv_t_jac[m][k] * inv_t_jac[j][r] +
                                    inv_t_jac[m][r] * inv_t_jac[j][k]);
                      }
                  grad_out[i][j] += tmp * inv_det;
                }
          }
        }
      return grad_out;
    }
  else
    {
      return BaseClass::get_gradient(q_point);
    }
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  get_divergence(const unsigned int q_point) const
{
#  ifdef DEBUG
  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  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 (this->data->element_type ==
      internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
    {
      if (!is_face &&
          this->cell_type == internal::MatrixFreeFunctions::cartesian)
        {
          // Cartesian cell
          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];
 
          // div * det(J^-1)
          divergence = this->gradients_quad[q_point] * inv_det;
          for (unsigned int d = 1; d < dim; ++d)
            divergence +=
              this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
        }
      else
        {
          // General cell
          // Derivatives are reordered for faces. Need to take this into account
          const VectorizedArrayType inv_det =
            determinant(
              this->jacobian[this->cell_type >
                                 internal::MatrixFreeFunctions::affine ?
                               q_point :
                               0]) *
            Number((is_face && dim == 2 && this->get_face_no() < 2) ? -1 : 1);
 
          // div * det(J^-1)
          divergence = this->gradients_quad[q_point] * inv_det;
          for (unsigned int d = 1; d < dim; ++d)
            divergence +=
              this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
        }
    }
  else
    {
      if (!is_face &&
          this->cell_type == internal::MatrixFreeFunctions::cartesian)
        {
          // Cartesian cell
          divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
          for (unsigned int d = 1; d < dim; ++d)
            divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
                          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];
          for (unsigned int e = 1; e < dim; ++e)
            divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
          for (unsigned int d = 1; d < dim; ++d)
            for (unsigned int e = 0; e < dim; ++e)
              divergence +=
                jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
        }
    }
  return divergence;
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, VectorizedArrayType>
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  get_symmetric_gradient(const unsigned int q_point) const
{
  // 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:
        Assert(false, ExcNotImplemented());
    }
  return SymmetricTensor<2, dim, VectorizedArrayType>(symmetrized);
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE
  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
  FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_curl(
    const unsigned int q_point) const
{
  // copy from generic function into dim-specialization function
  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
  switch (dim)
    {
      case 1:
        Assert(false,
               ExcMessage(
                 "Computing the curl in 1d is not a useful operation"));
        break;
      case 2:
        curl[0] = grad[1][0] - grad[0][1];
        break;
      case 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];
        break;
      default:
        Assert(false, ExcNotImplemented());
    }
  return curl;
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  get_hessian_diagonal(const unsigned int q_point) const
{
  return BaseClass::get_hessian_diagonal(q_point);
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<3, dim, VectorizedArrayType>
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_hessian(
  const unsigned int q_point) const
{
#  ifdef DEBUG
  Assert(this->hessians_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  return BaseClass::get_hessian(q_point);
}
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  submit_value(const Tensor<1, dim, VectorizedArrayType> val_in,
               const unsigned int                        q_point)
{
  if (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"));
#  ifdef DEBUG
      Assert(this->is_reinitialized, ExcNotInitialized());
      this->values_quad_submitted = true;
#  endif
 
      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)
            this->values_quad[comp * nqp + q_point] =
              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];
          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
          // 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)));
 
          // J^T * u * factor
          for (unsigned int comp = 0; comp < n_components; ++comp)
            {
              this->values_quad[comp * nqp + q_point] =
                val_in[0] * jac[0][comp] * fac;
              for (unsigned int e = 1; e < dim; ++e)
                this->values_quad[comp * nqp + q_point] +=
                  val_in[e] * jac[e][comp] * fac;
            }
        }
    }
  else
    {
      // No Piola transform
      BaseClass::submit_value(val_in, q_point);
    }
}
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  submit_gradient(const Tensor<2, dim, VectorizedArrayType> grad_in,
                  const unsigned int                        q_point)
{
  if (this->data->element_type ==
      internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
    {
      // Piola transform is required
 
#  ifdef DEBUG
      Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
      AssertIndexRange(q_point, this->n_quadrature_points);
      Assert(this->J_value != nullptr,
             internal::ExcMatrixFreeAccessToUninitializedMappingField(
               "update_gradients"));
      Assert(this->jacobian != nullptr,
             internal::ExcMatrixFreeAccessToUninitializedMappingField(
               "update_gradients"));
#  ifdef DEBUG
      this->gradients_quad_submitted = true;
#  endif
 
      const std::size_t nqp = this->n_quadrature_points;
      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)
              this->gradients_quad[(comp * dim + d) * nqp + q_point] =
                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
          for (unsigned int comp = 0; comp < n_components; ++comp)
            for (unsigned int d = 0; d < dim; ++d)
              {
                VectorizedArrayType tmp = 0;
                for (unsigned int f = 0; f < dim; ++f)
                  for (unsigned int e = 0; e < dim; ++e)
                    tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e];
 
                this->gradients_quad[(comp * dim + d) * nqp + q_point] =
                  tmp * 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];
          const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
 
          // 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));
 
          VectorizedArrayType tmp;
          // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
          for (unsigned int comp = 0; comp < n_components; ++comp)
            for (unsigned int d = 0; d < dim; ++d)
              {
                tmp = 0;
                for (unsigned int f = 0; f < dim; ++f)
                  for (unsigned int e = 0; e < dim; ++e)
                    tmp += t_jac[comp][f] * inv_t_jac[e][d] * grad_in[f][e];
 
                this->gradients_quad[(comp * dim + d) * nqp + q_point] =
                  tmp * fac;
              }
 
          // Contribution from values
          {
            // Diagonal part of jac_grad
 
            // Add jac_grad * J^{-1} * values * factor
            // -(J^{-T} * jac_grad * J^{-1} * J * values * factor)
            for (unsigned int f = 0; f < dim; ++f)
              {
                tmp = 0;
                for (unsigned int i = 0; i < dim; ++i)
                  for (unsigned int j = 0; j < dim; ++j)
                    {
                      tmp += inv_t_jac[j][f] * jac_grad[f][i] * grad_in[i][j];
                      for (unsigned int m = 0; m < dim; ++m)
                        for (unsigned int k = 0; k < dim; ++k)
                          tmp -= inv_t_jac[m][k] * jac_grad[k][m] *
                                 inv_t_jac[j][k] * t_jac[f][i] * grad_in[i][j];
                    }
                this->values_from_gradients_quad[f * nqp + q_point] = tmp * fac;
              }
          }
 
          {
            // Off-diagonal part of jac_grad
 
            // Add jac_grad * J^{-1} * values * factor
            for (unsigned int r = 0, f = dim; r < dim; ++r)
              for (unsigned int k = r + 1; k < dim; ++k, ++f)
                {
                  tmp = jac_grad[f][0] * inv_t_jac[0][k] * grad_in[0][0];
                  for (unsigned int j = 1; j < dim; ++j)
                    tmp += jac_grad[f][0] * inv_t_jac[j][k] * grad_in[0][j];
                  for (unsigned int i = 1; i < dim; ++i)
                    for (unsigned int j = 0; j < dim; ++j)
                      tmp += jac_grad[f][i] * inv_t_jac[j][k] * grad_in[i][j];
                  this->values_from_gradients_quad[r * nqp + q_point] +=
                    tmp * fac;
 
                  tmp = jac_grad[f][0] * inv_t_jac[0][r] * grad_in[0][0];
                  for (unsigned int j = 1; j < dim; ++j)
                    tmp += jac_grad[f][0] * inv_t_jac[j][r] * grad_in[0][j];
                  for (unsigned int i = 1; i < dim; ++i)
                    for (unsigned int j = 0; j < dim; ++j)
                      tmp += jac_grad[f][i] * inv_t_jac[j][r] * grad_in[i][j];
                  this->values_from_gradients_quad[k * nqp + q_point] +=
                    tmp * fac;
                }
 
            // -(J^{-T} * jac_grad * J^{-1} * J * values * factor)
            for (unsigned int n = 0; n < dim; ++n)
              {
                tmp = 0;
                for (unsigned int r = 0, f = dim; r < dim; ++r)
                  for (unsigned int k = r + 1; k < dim; ++k, ++f)
                    for (unsigned int i = 0; i < dim; ++i)
                      for (unsigned int j = 0; j < dim; ++j)
                        for (unsigned int m = 0; m < dim; ++m)
                          tmp += jac_grad[f][m] * t_jac[n][i] * grad_in[i][j] *
                                 (inv_t_jac[m][k] * inv_t_jac[j][r] +
                                  inv_t_jac[m][r] * inv_t_jac[j][k]);
 
                this->values_from_gradients_quad[n * nqp + q_point] -=
                  tmp * fac;
              }
          }
        }
    }
  else
    {
      BaseClass::submit_gradient(grad_in, q_point);
    }
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  submit_gradient(
    const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
    const unsigned int                                        q_point)
{
  if (this->data->element_type ==
      internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
    {
      // Piola transform is required
      const Tensor<2, dim, VectorizedArrayType> &grad = grad_in;
      FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
        submit_gradient(grad, q_point);
    }
  else
    {
      BaseClass::submit_gradient(grad_in, q_point);
    }
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  submit_divergence(const VectorizedArrayType div_in,
                    const unsigned int        q_point)
{
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  const std::size_t nqp = this->n_quadrature_points;
  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)
        {
          this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
          for (unsigned int e = d + 1; e < dim; ++e)
            {
              this->gradients_quad[(dim * d + e) * nqp + q_point] =
                VectorizedArrayType();
              this->gradients_quad[(dim * e + d) * nqp + q_point] =
                VectorizedArrayType();
            }
        }
      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)
            {
              this->gradients_quad[(d * dim + d) * nqp + q_point] =
                (fac * this->jacobian[0][d][d]);
              for (unsigned int e = d + 1; e < dim; ++e)
                {
                  this->gradients_quad[(d * dim + e) * nqp + q_point] =
                    VectorizedArrayType();
                  this->gradients_quad[(e * dim + d) * nqp + q_point] =
                    VectorizedArrayType();
                }
            }
        }
      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)
                this->gradients_quad[(d * dim + e) * nqp + q_point] =
                  jac[d][e] * fac;
            }
        }
    }
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  submit_symmetric_gradient(
    const SymmetricTensor<2, dim, VectorizedArrayType> sym_grad,
    const unsigned int                                 q_point)
{
  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
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  Assert(this->jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  const std::size_t nqp = this->n_quadrature_points;
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
      for (unsigned int d = 0; d < dim; ++d)
        this->gradients_quad[(d * dim + d) * nqp + q_point] =
          (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][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;
            this->gradients_quad[(e * dim + d) * nqp + q_point] =
              value * this->jacobian[0][d][d];
            this->gradients_quad[(d * dim + e) * nqp + q_point] =
              value * this->jacobian[0][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];
            this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
          }
    }
}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::submit_curl(
  const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl,
  const unsigned int                                       q_point)
{
  Tensor<2, dim, VectorizedArrayType> grad;
  switch (dim)
    {
      case 1:
        Assert(false,
               ExcMessage(
                 "Testing by the curl in 1d is not a useful operation"));
        break;
      case 2:
        grad[1][0] = curl[0];
        grad[0][1] = -curl[0];
        break;
      case 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];
        break;
      default:
        Assert(false, ExcNotImplemented());
    }
  submit_gradient(grad, q_point);
}
 
 
/*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
  FEEvaluationAccess(
    const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
    const unsigned int                                dof_no,
    const unsigned int                                first_selected_component,
    const unsigned int                                quad_no,
    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<1, 1, Number, is_face, VectorizedArrayType>(
      matrix_free,
      dof_no,
      first_selected_component,
      quad_no,
      fe_degree,
      n_q_points,
      is_interior_face,
      active_fe_index,
      active_quad_index,
      face_type)
{}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
  FEEvaluationAccess(
    const Mapping<1> &      mapping,
    const FiniteElement<1> &fe,
    const Quadrature<1> &   quadrature,
    const UpdateFlags       update_flags,
    const unsigned int      first_selected_component,
    const FEEvaluationData<1, VectorizedArrayType, is_face> *other)
  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
      mapping,
      fe,
      quadrature,
      update_flags,
      first_selected_component,
      other)
{}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
  FEEvaluationAccess(
    const FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType> &other)
  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(other)
{}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType> &
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::operator=(
  const FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType> &other)
{
  this->FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>::operator=(
    other);
  return *this;
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::get_dof_value(
  const unsigned int dof) const
{
  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
  return this->values_dofs[dof];
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::get_value(
  const unsigned int q_point) const
{
#  ifdef DEBUG
  Assert(this->values_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
  return this->values_quad[q_point];
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<1, 1, VectorizedArrayType>
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::get_gradient(
  const unsigned int q_point) const
{
  // could use the base class gradient, but that involves too many inefficient
  // initialization operations on tensors
 
#  ifdef DEBUG
  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
 
  const Tensor<2, 1, VectorizedArrayType> &jac =
    this->cell_type == internal::MatrixFreeFunctions::general ?
      this->jacobian[q_point] :
      this->jacobian[0];
 
  Tensor<1, 1, VectorizedArrayType> grad_out;
  grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
 
  return grad_out;
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::get_divergence(
  const unsigned int q_point) const
{
  return get_gradient(q_point)[0];
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
  get_normal_derivative(const unsigned int q_point) const
{
  return BaseClass::get_normal_derivative(q_point)[0];
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<2, 1, VectorizedArrayType>
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::get_hessian(
  const unsigned int q_point) const
{
  return BaseClass::get_hessian(q_point)[0];
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<1, 1, VectorizedArrayType>
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
  get_hessian_diagonal(const unsigned int q_point) const
{
  return BaseClass::get_hessian_diagonal(q_point)[0];
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::get_laplacian(
  const unsigned int q_point) const
{
  return BaseClass::get_laplacian(q_point)[0];
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void DEAL_II_ALWAYS_INLINE
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
{
#  ifdef DEBUG
  this->dof_values_initialized = true;
  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
#  endif
  this->values_dofs[dof] = val_in;
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::submit_value(
  const VectorizedArrayType val_in,
  const unsigned int        q_point)
{
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
#  ifdef DEBUG
  this->values_quad_submitted = true;
#  endif
 
  if (this->cell_type == internal::MatrixFreeFunctions::general)
    {
      const VectorizedArrayType JxW = this->J_value[q_point];
      this->values_quad[q_point]    = val_in * JxW;
    }
  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_point];
      this->values_quad[q_point] = val_in * JxW;
    }
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::submit_value(
  const Tensor<1, 1, VectorizedArrayType> val_in,
  const unsigned int                      q_point)
{
  submit_value(val_in[0], q_point);
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::submit_gradient(
  const Tensor<1, 1, VectorizedArrayType> grad_in,
  const unsigned int                      q_point)
{
  submit_gradient(grad_in[0], q_point);
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::submit_gradient(
  const VectorizedArrayType grad_in,
  const unsigned int        q_point)
{
#  ifdef DEBUG
  Assert(this->is_reinitialized, ExcNotInitialized());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  const Tensor<2, 1, VectorizedArrayType> jac =
    this->cell_type == internal::MatrixFreeFunctions::general ?
      this->jacobian[q_point] :
      this->jacobian[0];
  const VectorizedArrayType JxW =
    this->cell_type == internal::MatrixFreeFunctions::general ?
      this->J_value[q_point] :
      this->J_value[0] * this->quadrature_weights[q_point];
 
  this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::submit_gradient(
  const Tensor<2, 1, VectorizedArrayType> grad_in,
  const unsigned int                      q_point)
{
  submit_gradient(grad_in[0][0], q_point);
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
  submit_normal_derivative(const VectorizedArrayType grad_in,
                           const unsigned int        q_point)
{
  Tensor<1, 1, VectorizedArrayType> grad;
  grad[0] = grad_in;
  BaseClass::submit_normal_derivative(grad, q_point);
}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
  submit_normal_derivative(const Tensor<1, 1, VectorizedArrayType> grad_in,
                           const unsigned int                      q_point)
{
  BaseClass::submit_normal_derivative(grad_in, q_point);
}
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::submit_hessian(
  const Tensor<2, 1, VectorizedArrayType> hessian_in,
  const unsigned int                      q_point)
{
  Tensor<1, 1, Tensor<2, 1, VectorizedArrayType>> hessian;
  hessian[0] = hessian_in;
  BaseClass::submit_hessian(hessian, q_point);
}
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
inline VectorizedArrayType
FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
  integrate_value() const
{
  return BaseClass::integrate_value()[0];
}
 
 
 
/*-------------------------- 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                                  quad_no,
               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,
              quad_no,
              fe_degree,
              static_n_q_points,
              true /*note: this is not a face*/,
              active_fe_index,
              active_quad_index)
  , 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_no,
               const unsigned int                                  quad_no,
               const unsigned int first_selected_component)
  : FEEvaluation(matrix_free,
                 dof_no,
                 quad_no,
                 first_selected_component,
                 matrix_free.get_cell_active_fe_index(range))
{}
 
 
 
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_no,
                           const unsigned int first_selected_component)
{
  (void)dof_no;
  (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."));
 
#  ifdef DEBUG
  // 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_no) + ", ";
          message += Utilities::int_to_string(this->quad_no) + ", ";
          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_no != numbers::invalid_unsigned_int)
        {
          if (static_cast<unsigned int>(fe_degree) ==
              this->data->data.front().fe_degree)
            {
              proposed_dof_comp = dof_no;
              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->quad_no;
          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_no != 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_no)
            correct_pos = " ^ ";
          else
            correct_pos = "   ";
          if (proposed_quad_comp != this->quad_no)
            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_no != numbers::invalid_unsigned_int)
        {
          message += ", " + Utilities::int_to_string(dof_no) + ", ";
          message += Utilities::int_to_string(this->quad_no);
          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_no != numbers::invalid_unsigned_int)
    AssertDimension(
      n_q_points,
      this->mapping_data->descriptor[this->active_quad_index].n_q_points);
#  endif
}
 
 
 
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) ==
      VectorizedArrayType::size())
    {
      DEAL_II_OPENMP_SIMD_PRAGMA
      for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
        this->cell_ids[i] = cell_index * VectorizedArrayType::size() + 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 * VectorizedArrayType::size() + i;
      for (; i < VectorizedArrayType::size(); ++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]];
 
#  ifdef DEBUG
  this->is_reinitialized           = true;
  this->dof_values_initialized     = false;
  this->values_quad_initialized    = false;
  this->gradients_quad_initialized = false;
  this->hessians_quad_initialized  = false;
#  endif
}
 
 
 
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, VectorizedArrayType::size()> &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 < VectorizedArrayType::size(); ++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 / VectorizedArrayType::size()));
    }
 
  // 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->mapping_data->jacobians[0].size() > 0)
        this_jacobian_data.resize_fast(2);
 
      if (this->mapping_data->JxW_values.size() > 0)
        this_J_value_data.resize_fast(1);
 
      if (this->mapping_data->jacobian_gradients[0].size() > 0)
        this_jacobian_gradients_data.resize_fast(1);
 
      if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
        this_jacobian_gradients_non_inverse_data.resize_fast(1);
 
      if (this->mapping_data->quadrature_points.size() > 0)
        this_quadrature_points_data.resize_fast(1);
    }
  else
    {
      if (this->mapping_data->jacobians[0].size() > 0)
        this_jacobian_data.resize_fast(this->n_quadrature_points);
 
      if (this->mapping_data->JxW_values.size() > 0)
        this_J_value_data.resize_fast(this->n_quadrature_points);
 
      if (this->mapping_data->jacobian_gradients[0].size() > 0)
        this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
 
      if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
        this_jacobian_gradients_non_inverse_data.resize_fast(
          this->n_quadrature_points);
 
      if (this->mapping_data->quadrature_points.size() > 0)
        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 < VectorizedArrayType::size(); ++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 / VectorizedArrayType::size();
      const unsigned int offsets =
        this->mapping_data->data_index_offsets[cell_batch_index];
      const unsigned int lane = cell_index % VectorizedArrayType::size();
 
      if (this->cell_type <=
          internal::MatrixFreeFunctions::GeometryType::affine)
        {
          // case that all cells are Cartesian or affine
          const unsigned int q = 0;
 
          if (this->mapping_data->JxW_values.size() > 0)
            this_J_value_data[q][v] =
              this->mapping_data->JxW_values[offsets + q][lane];
 
          if (this->mapping_data->jacobians[0].size() > 0)
            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];
 
          if (this->mapping_data->jacobian_gradients[0].size() > 0)
            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];
 
          if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
            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 (this->mapping_data->quadrature_points.size() > 0)
            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;
 
              if (this->mapping_data->JxW_values.size() > 0)
                this_J_value_data[q][v] =
                  this->mapping_data->JxW_values[offsets + q_src][lane];
 
              if (this->mapping_data->jacobians[0].size() > 0)
                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];
 
              if (this->mapping_data->jacobian_gradients[0].size() > 0)
                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];
 
              if (this->mapping_data->jacobian_gradients_non_inverse[0].size() >
                  0)
                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 (this->mapping_data->quadrature_points.size() > 0)
                {
                  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];
                    }
                }
            }
        }
    }
 
#  ifdef DEBUG
  this->is_reinitialized           = true;
  this->dof_values_initialized     = false;
  this->values_quad_initialized    = false;
  this->gradients_quad_initialized = false;
  this->hessians_quad_initialized  = false;
#  endif
}
 
 
 
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);
 
#  ifdef DEBUG
  this->is_reinitialized = true;
#  endif
}
 
 
 
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);
 
#  ifdef DEBUG
  this->is_reinitialized = true;
#  endif
}
 
 
 
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 bool evaluate_values,
                                            const bool evaluate_gradients,
                                            const bool evaluate_hessians)
{
#  ifdef DEBUG
  Assert(this->dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  evaluate(this->values_dofs,
           evaluate_values,
           evaluate_gradients,
           evaluate_hessians);
}
 
 
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)
{
#  ifdef DEBUG
  Assert(this->dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  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 bool evaluate_values,
                                            const bool evaluate_gradients,
                                            const bool evaluate_hessians)
{
  const EvaluationFlags::EvaluationFlags flag =
    ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
    ((evaluate_gradients) ? EvaluationFlags::gradients :
                            EvaluationFlags::nothing) |
    ((evaluate_hessians) ? EvaluationFlags::hessians :
                           EvaluationFlags::nothing);
 
  evaluate(values_array, 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>::
  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 (fe_degree > -1)
    {
      SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
        evaluate(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);
    }
 
#  ifdef DEBUG
  if (evaluation_flag_actual & EvaluationFlags::values)
    this->values_quad_initialized = true;
  if (evaluation_flag_actual & EvaluationFlags::gradients)
    this->gradients_quad_initialized = true;
  if (evaluation_flag_actual & EvaluationFlags::hessians)
    this->hessians_quad_initialized = true;
#  endif
}
 
 
 
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 bool        evaluate_values,
                                        const bool        evaluate_gradients,
                                        const bool        evaluate_hessians)
{
  const EvaluationFlags::EvaluationFlags flag =
    ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
    ((evaluate_gradients) ? EvaluationFlags::gradients :
                            EvaluationFlags::nothing) |
    ((evaluate_hessians) ? EvaluationFlags::hessians :
                           EvaluationFlags::nothing);
 
  gather_evaluate(input_vector, flag);
}
 
 
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<typename VectorType::value_type, Number>::value &&
        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 bool integrate_values,
                                             const bool integrate_gradients)
{
  integrate(integrate_values, integrate_gradients, this->values_dofs);
 
#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif
}
 
 
 
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);
 
#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif
}
 
 
 
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 bool integrate_values,
                                             const bool integrate_gradients,
                                             VectorizedArrayType *values_array)
{
  EvaluationFlags::EvaluationFlags flag =
    (integrate_values ? EvaluationFlags::values : EvaluationFlags::nothing) |
    (integrate_gradients ? EvaluationFlags::gradients :
                           EvaluationFlags::nothing);
  integrate(flag, values_array);
}
 
 
 
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)
{
#  ifdef DEBUG
  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());
#  endif
  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 (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 (fe_degree > -1)
    {
      SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
        integrate(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);
    }
 
#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif
}
 
 
 
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 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
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 {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_no,
    const unsigned int                                  quad_no,
    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_no,
              first_selected_component,
              quad_no,
              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_no,
    const unsigned int                                  quad_no,
    const unsigned int first_selected_component)
  : FEFaceEvaluation(matrix_free,
                     is_interior_face,
                     dof_no,
                     quad_no,
                     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 * VectorizedArrayType::size() + i;
  for (; i < VectorizedArrayType::size(); ++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];
    }
 
#  ifdef DEBUG
  this->is_reinitialized           = true;
  this->dof_values_initialized     = false;
  this->values_quad_initialized    = false;
  this->gradients_quad_initialized = false;
  this->hessians_quad_initialized  = false;
#  endif
}
 
 
 
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->quad_no <
      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;
 
  constexpr unsigned int n_lanes = VectorizedArrayType::size();
 
  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] = static_cast<std::uint8_t>(-1);
              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;
          unsigned int face_orientation          = faces.face_orientation % 8;
          if (face_identifies_as_interior != orientation_interior_face)
            {
              constexpr std::array<std::uint8_t, 8> table{
                {0, 1, 2, 3, 6, 5, 4, 7}};
              face_orientation = table[face_orientation];
            }
          this->face_orientations[i] = face_orientation;
        }
    }
  else
    {
      this->face_orientations[0] = 0;
      this->face_numbers[0]      = face_number;
      for (unsigned int i = 0; i < n_lanes; ++i)
        this->cell_ids[i] = cell_index * n_lanes + i;
      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->quad_no]
      .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->quad_no]
                     .JxW_values.size());
  this->J_value = &this->matrix_free->get_mapping_info()
                     .face_data_by_cells[this->quad_no]
                     .JxW_values[offsets];
  this->normal_vectors = &this->matrix_free->get_mapping_info()
                            .face_data_by_cells[this->quad_no]
                            .normal_vectors[offsets];
  this->jacobian = &this->matrix_free->get_mapping_info()
                      .face_data_by_cells[this->quad_no]
                      .jacobians[!this->is_interior_face()][offsets];
  this->normal_x_jacobian =
    &this->matrix_free->get_mapping_info()
       .face_data_by_cells[this->quad_no]
       .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->quad_no]
        .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->quad_no]
                         .quadrature_point_offsets.size());
      this->quadrature_points = this->matrix_free->get_mapping_info()
                                  .face_data_by_cells[this->quad_no]
                                  .quadrature_points.data() +
                                this->matrix_free->get_mapping_info()
                                  .face_data_by_cells[this->quad_no]
                                  .quadrature_point_offsets[index];
    }
 
#  ifdef DEBUG
  this->is_reinitialized           = true;
  this->dof_values_initialized     = false;
  this->values_quad_initialized    = false;
  this->gradients_quad_initialized = false;
  this->hessians_quad_initialized  = false;
#  endif
}
 
 
 
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 bool evaluate_values,
                                                const bool evaluate_gradients)
{
#  ifdef DEBUG
  Assert(this->dof_values_initialized, ExcNotInitialized());
#  endif
 
  evaluate(this->values_dofs, evaluate_values, evaluate_gradients);
}
 
 
 
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)
{
#  ifdef DEBUG
  Assert(this->dof_values_initialized, ExcNotInitialized());
#  endif
 
  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 bool evaluate_values,
                                                const bool evaluate_gradients)
{
  const EvaluationFlags::EvaluationFlags flag =
    ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
    ((evaluate_gradients) ? EvaluationFlags::gradients :
                            EvaluationFlags::nothing);
 
  evaluate(values_array, 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 (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);
 
#  ifdef DEBUG
  if (evaluation_flag_actual & EvaluationFlags::values)
    this->values_quad_initialized = true;
  if (evaluation_flag_actual & EvaluationFlags::gradients)
    this->gradients_quad_initialized = true;
  if ((evaluation_flag_actual & EvaluationFlags::hessians) != 0u)
    this->hessians_quad_initialized = true;
#  endif
}
 
 
 
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)
{
  integrate(integration_flag, this->values_dofs);
 
#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif
}
 
 
 
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 bool integrate_values,
                                                 const bool integrate_gradients)
{
  integrate(integrate_values, integrate_gradients, this->values_dofs);
 
#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif
}
 
 
 
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 bool integrate_values,
                                                 const bool integrate_gradients,
                                                 VectorizedArrayType
                                                   *values_array)
{
  const EvaluationFlags::EvaluationFlags flag =
    ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
    ((integrate_gradients) ? EvaluationFlags::gradients :
                             EvaluationFlags::nothing);
 
  integrate(flag, values_array);
}
 
 
 
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)
{
  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 (fe_degree > -1)
    internal::FEFaceEvaluationImplIntegrateSelector<dim, VectorizedArrayType>::
      template run<fe_degree, n_q_points_1d>(n_components,
                                             integration_flag_actual,
                                             values_array,
                                             *this);
  else
    internal::FEFaceEvaluationFactory<dim, VectorizedArrayType>::integrate(
      n_components, integration_flag_actual, values_array, *this);
}
 
 
 
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 bool        evaluate_values,
                                        const bool        evaluate_gradients)
{
  const EvaluationFlags::EvaluationFlags flag =
    ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
    ((evaluate_gradients) ? EvaluationFlags::gradients :
                            EvaluationFlags::nothing);
 
  gather_evaluate(input_vector, flag);
}
 
 
 
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 (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);
    }
 
#  ifdef DEBUG
  if (evaluation_flag & EvaluationFlags::values)
    this->values_quad_initialized = true;
  if (evaluation_flag & EvaluationFlags::gradients)
    this->gradients_quad_initialized = true;
  if (evaluation_flag & EvaluationFlags::hessians)
    this->hessians_quad_initialized = true;
#  endif
}
 
 
 
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 (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 {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;
}
 
 
 
/*------------------------- end FEFaceEvaluation ------------------------- */
 
 
#endif // ifndef DOXYGEN
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif