fe_evaluation.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2020 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_kernels.h>
#include <deal.II/matrix_free/evaluation_selector.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>
 
 
DEAL_II_NAMESPACE_OPEN
 
 
 
namespace internal
{
  DeclException0(ExcAccessToUninitializedField);
}
 
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 FEEvaluation;
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face                 = false,
          typename VectorizedArrayType = VectorizedArray<Number>>
class FEEvaluationBase
{
  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_;
 
  ~FEEvaluationBase();
 
  unsigned int
  get_mapping_data_index_offset() const;
 
  internal::MatrixFreeFunctions::GeometryType
  get_cell_type() const;
 
  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
  get_shape_info() const;
 
 
  template <typename VectorType>
  void
  read_dof_values(const VectorType &src, const unsigned int first_index = 0);
 
  template <typename VectorType>
  void
  read_dof_values_plain(const VectorType & src,
                        const unsigned int first_index = 0);
 
  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;
 
 
  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);
 
  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;
 
  VectorizedArrayType
  JxW(const unsigned int q_index) const;
 
  DEAL_II_DEPRECATED void
  fill_JxW_values(AlignedVector<VectorizedArrayType> &JxW_values) const;
 
  Tensor<2, dim, VectorizedArrayType>
  inverse_jacobian(const unsigned int q_index) const;
 
  Tensor<1, dim, VectorizedArrayType>
  get_normal_vector(const unsigned int q_point) const;
 
  VectorizedArrayType
  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const;
 
  template <typename T>
  std::array<T, VectorizedArrayType::size()>
  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
                   &array) const;
 
  std::array<unsigned int, VectorizedArrayType::size()>
  get_cell_ids() const;
 
 
  const VectorizedArrayType *
  begin_dof_values() const;
 
  VectorizedArrayType *
  begin_dof_values();
 
  const VectorizedArrayType *
  begin_values() const;
 
  VectorizedArrayType *
  begin_values();
 
  const VectorizedArrayType *
  begin_gradients() const;
 
  VectorizedArrayType *
  begin_gradients();
 
  const VectorizedArrayType *
  begin_hessians() const;
 
  VectorizedArrayType *
  begin_hessians();
 
  const std::vector<unsigned int> &
  get_internal_dof_numbering() const;
 
  ArrayView<VectorizedArrayType>
  get_scratch_data() 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);
 
  template <int n_components_other>
  FEEvaluationBase(const Mapping<dim> &      mapping,
                   const FiniteElement<dim> &fe,
                   const Quadrature<1> &     quadrature,
                   const UpdateFlags         update_flags,
                   const unsigned int        first_selected_component,
                   const FEEvaluationBase<dim,
                                          n_components_other,
                                          Number,
                                          is_face,
                                          VectorizedArrayType> *other);
 
  FEEvaluationBase(const FEEvaluationBase &other);
 
  FEEvaluationBase &
  operator=(const FEEvaluationBase &other);
 
  template <typename VectorType, typename VectorOperation>
  void
  read_write_operation(const VectorOperation &operation,
                       VectorType *           vectors[],
                       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,
    VectorType *                                    vectors[],
    const std::bitset<VectorizedArrayType::size()> &mask) const;
 
  template <typename VectorType, typename VectorOperation>
  void
  read_write_operation_global(const VectorOperation &operation,
                              VectorType *           vectors[]) const;
 
  AlignedVector<VectorizedArrayType> *scratch_data_array;
 
  VectorizedArrayType *scratch_data;
 
  VectorizedArrayType *values_dofs[n_components];
 
  VectorizedArrayType *values_quad[n_components];
 
  VectorizedArrayType *gradients_quad[n_components][dim];
 
  VectorizedArrayType *hessians_quad[n_components][(dim * (dim + 1)) / 2];
 
  const unsigned int quad_no;
 
  const unsigned int n_fe_components;
 
  const unsigned int active_fe_index;
 
  const unsigned int active_quad_index;
 
  const unsigned int n_quadrature_points;
 
  const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info;
 
  const internal::MatrixFreeFunctions::DoFInfo *dof_info;
 
  const internal::MatrixFreeFunctions::MappingInfoStorage<
    (is_face ? dim - 1 : dim),
    dim,
    Number,
    VectorizedArrayType> *mapping_data;
 
  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *data;
 
  const Tensor<2, dim, VectorizedArrayType> *jacobian;
 
  const VectorizedArrayType *J_value;
 
  const Tensor<1, dim, VectorizedArrayType> *normal_vectors;
 
  const Tensor<1, dim, VectorizedArrayType> *normal_x_jacobian;
 
  const Number *quadrature_weights;
 
  unsigned int cell;
 
  bool is_interior_face;
 
  internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index;
 
  unsigned int face_no;
 
  unsigned int face_orientation;
 
  unsigned int subface_index;
 
  internal::MatrixFreeFunctions::GeometryType cell_type;
 
  bool dof_values_initialized;
 
  bool values_quad_initialized;
 
  bool gradients_quad_initialized;
 
  bool hessians_quad_initialized;
 
  bool values_quad_submitted;
 
  bool gradients_quad_submitted;
 
  std::shared_ptr<internal::MatrixFreeFunctions::
                    MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
    mapped_geometry;
 
  const unsigned int first_selected_component;
 
  mutable std::vector<types::global_dof_index> local_dof_indices;
 
private:
  void
  set_data_pointers();
 
  // Make other FEEvaluationBase as well as FEEvaluation objects friends.
  template <int, int, typename, bool, typename>
  friend class FEEvaluationBase;
  template <int, int, int, int, typename, typename>
  friend class FEEvaluation;
};
 
 
 
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);
 
  template <int n_components_other>
  FEEvaluationAccess(const Mapping<dim> &      mapping,
                     const FiniteElement<dim> &fe,
                     const Quadrature<1> &     quadrature,
                     const UpdateFlags         update_flags,
                     const unsigned int        first_selected_component,
                     const FEEvaluationBase<dim,
                                            n_components_other,
                                            Number,
                                            is_face,
                                            VectorizedArrayType> *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>;
  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);
 
  Tensor<2, dim, VectorizedArrayType>
  get_hessian(unsigned int q_point) const;
 
  gradient_type
  get_hessian_diagonal(const unsigned int q_point) const;
 
  value_type
  get_laplacian(const unsigned int q_point) const;
 
  value_type
  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);
 
  template <int n_components_other>
  FEEvaluationAccess(const Mapping<dim> &      mapping,
                     const FiniteElement<dim> &fe,
                     const Quadrature<1> &     quadrature,
                     const UpdateFlags         update_flags,
                     const unsigned int        first_selected_component,
                     const FEEvaluationBase<dim,
                                            n_components_other,
                                            Number,
                                            is_face,
                                            VectorizedArrayType> *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>;
 
  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_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);
 
  template <int n_components_other>
  FEEvaluationAccess(const Mapping<dim> &      mapping,
                     const FiniteElement<dim> &fe,
                     const Quadrature<1> &     quadrature,
                     const UpdateFlags         update_flags,
                     const unsigned int        first_selected_component,
                     const FEEvaluationBase<dim,
                                            n_components_other,
                                            Number,
                                            is_face,
                                            VectorizedArrayType> *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>;
  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_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_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);
 
  Tensor<2, 1, VectorizedArrayType>
  get_hessian(unsigned int q_point) const;
 
  gradient_type
  get_hessian_diagonal(const unsigned int q_point) const;
 
  value_type
  get_laplacian(const unsigned int q_point) const;
 
  value_type
  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);
 
  template <int n_components_other>
  FEEvaluationAccess(const Mapping<1> &      mapping,
                     const FiniteElement<1> &fe,
                     const Quadrature<1> &   quadrature,
                     const UpdateFlags       update_flags,
                     const unsigned int      first_selected_component,
                     const FEEvaluationBase<1,
                                            n_components_other,
                                            Number,
                                            is_face,
                                            VectorizedArrayType> *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);
 
  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);
 
  template <int n_components_other>
  FEEvaluation(const FiniteElement<dim> &                   fe,
               const FEEvaluationBase<dim,
                                      n_components_other,
                                      Number,
                                      false,
                                      VectorizedArrayType> &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);
 
  template <typename DoFHandlerType, bool level_dof_access>
  void
  reinit(const TriaIterator<DoFCellAccessor<DoFHandlerType, level_dof_access>>
           &cell);
 
  void
  reinit(const typename Triangulation<dim>::cell_iterator &cell);
 
  void
  evaluate(const bool evaluate_values,
           const bool evaluate_gradients,
           const bool evaluate_hessians = false);
 
  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 bool        evaluate_values,
                  const bool        evaluate_gradients,
                  const bool        evaluate_hessians = false);
 
  void
  integrate(const bool integrate_values, const bool integrate_gradients);
 
  void
  integrate(const bool           integrate_values,
            const bool           integrate_gradients,
            VectorizedArrayType *values_array);
 
  template <typename VectorType>
  void
  integrate_scatter(const bool  integrate_values,
                    const bool  integrate_gradients,
                    VectorType &output_vector);
 
  Point<dim, VectorizedArrayType>
  quadrature_point(const unsigned int q_point) 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);
 
  void
  reinit(const unsigned int face_batch_number);
 
  void
  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
 
  void
  evaluate(const bool evaluate_values, const bool evaluate_gradients);
 
  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 bool        evaluate_values,
                  const bool        evaluate_gradients);
 
  void
  integrate(const bool integrate_values, const bool integrate_gradients);
 
  void
  integrate(const bool           integrate_values,
            const bool           integrate_gradients,
            VectorizedArrayType *values_array);
 
  template <typename VectorType>
  void
  integrate_scatter(const bool  integrate_values,
                    const bool  integrate_gradients,
                    VectorType &output_vector);
 
  Point<dim, VectorizedArrayType>
  quadrature_point(const unsigned int q_point) 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
 
 
 
/*----------------------- 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> &data_in,
                   const unsigned int                                  dof_no,
                   const unsigned int first_selected_component,
                   const unsigned int quad_no_in,
                   const unsigned int fe_degree,
                   const unsigned int n_q_points,
                   const bool         is_interior_face)
  : scratch_data_array(data_in.acquire_scratch_data())
  , quad_no(quad_no_in)
  , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
  , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
                      data_in.get_dof_info(dof_no).fe_index_from_degree(
                        first_selected_component,
                        fe_degree) :
                      0)
  , active_quad_index(fe_degree != numbers::invalid_unsigned_int ?
                        (is_face ? data_in.get_mapping_info()
                                     .face_data[quad_no_in]
                                     .quad_index_from_n_q_points(n_q_points) :
                                   data_in.get_mapping_info()
                                     .cell_data[quad_no_in]
                                     .quad_index_from_n_q_points(n_q_points)) :
                        0)
  , n_quadrature_points(fe_degree != numbers::invalid_unsigned_int ?
                          n_q_points :
                          (is_face ? data_in
                                       .get_shape_info(dof_no,
                                                       quad_no_in,
                                                       active_fe_index,
                                                       active_quad_index)
                                       .n_q_points_face :
                                     data_in
                                       .get_shape_info(dof_no,
                                                       quad_no_in,
                                                       active_fe_index,
                                                       active_quad_index)
                                       .n_q_points))
  , matrix_info(&data_in)
  , dof_info(&data_in.get_dof_info(dof_no))
  , mapping_data(
      internal::MatrixFreeFunctions::
        MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
          data_in.get_mapping_info(),
          quad_no))
  , data(&data_in.get_shape_info(
      dof_no,
      quad_no_in,
      dof_info->component_to_base_index[first_selected_component],
      active_fe_index,
      active_quad_index))
  , jacobian(nullptr)
  , J_value(nullptr)
  , normal_vectors(nullptr)
  , normal_x_jacobian(nullptr)
  , quadrature_weights(
      mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
  , cell(numbers::invalid_unsigned_int)
  , is_interior_face(is_interior_face)
  , dof_access_index(
      is_face ?
        (is_interior_face ?
           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
  , cell_type(internal::MatrixFreeFunctions::general)
  , dof_values_initialized(false)
  , values_quad_initialized(false)
  , gradients_quad_initialized(false)
  , hessians_quad_initialized(false)
  , values_quad_submitted(false)
  , gradients_quad_submitted(false)
  , first_selected_component(first_selected_component)
{
  set_data_pointers();
  Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
  AssertDimension(matrix_info->get_task_info().vectorization_length,
                  VectorizedArrayType::size());
  AssertDimension((is_face ? data->n_q_points_face : data->n_q_points),
                  n_quadrature_points);
  AssertDimension(n_quadrature_points,
                  mapping_data->descriptor[active_quad_index].n_q_points);
  Assert(
    dof_info->start_components.back() == 1 ||
      static_cast<int>(n_components_) <=
        static_cast<int>(
          dof_info->start_components
            [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(
        dof_info->start_components
          [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 -
        dof_info->start_components
          [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>
template <int n_components_other>
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 FEEvaluationBase<dim,
                                          n_components_other,
                                          Number,
                                          is_face,
                                          VectorizedArrayType> *other)
  : scratch_data_array(new AlignedVector<VectorizedArrayType>())
  , quad_no(numbers::invalid_unsigned_int)
  , n_fe_components(n_components_)
  , active_fe_index(numbers::invalid_unsigned_int)
  , active_quad_index(numbers::invalid_unsigned_int)
  , n_quadrature_points(
      Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
  , matrix_info(nullptr)
  , dof_info(nullptr)
  , mapping_data(nullptr)
  ,
  // select the correct base element from the given FE component
  data(new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
    quadrature,
    fe,
    fe.component_to_base_index(first_selected_component).first))
  , jacobian(nullptr)
  , J_value(nullptr)
  , normal_vectors(nullptr)
  , normal_x_jacobian(nullptr)
  , quadrature_weights(nullptr)
  , cell(0)
  , cell_type(internal::MatrixFreeFunctions::general)
  , is_interior_face(true)
  , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
  , dof_values_initialized(false)
  , values_quad_initialized(false)
  , gradients_quad_initialized(false)
  , hessians_quad_initialized(false)
  , values_quad_submitted(false)
  , gradients_quad_submitted(false)
  ,
  // keep the number of the selected component within the current base element
  // for reading dof values
  first_selected_component(first_selected_component)
{
  set_data_pointers();
 
  Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
         ExcInternalError());
  if (other != nullptr &&
      other->mapped_geometry->get_quadrature() == quadrature)
    mapped_geometry = other->mapped_geometry;
  else
    mapped_geometry =
      std::make_shared<internal::MatrixFreeFunctions::
                         MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
        mapping, quadrature, update_flags);
  cell = 0;
 
  mapping_data = &mapped_geometry->get_data_storage();
  jacobian     = mapped_geometry->get_data_storage().jacobians[0].begin();
  J_value      = mapped_geometry->get_data_storage().JxW_values.begin();
 
  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;
}
 
 
 
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)
  : scratch_data_array(other.matrix_info == nullptr ?
                         new AlignedVector<VectorizedArrayType>() :
                         other.matrix_info->acquire_scratch_data())
  , quad_no(other.quad_no)
  , n_fe_components(other.n_fe_components)
  , active_fe_index(other.active_fe_index)
  , active_quad_index(other.active_quad_index)
  , n_quadrature_points(other.n_quadrature_points)
  , matrix_info(other.matrix_info)
  , dof_info(other.dof_info)
  , mapping_data(other.mapping_data)
  , data(other.matrix_info == nullptr ?
           new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
             *other.data) :
           other.data)
  , jacobian(nullptr)
  , J_value(nullptr)
  , normal_vectors(nullptr)
  , normal_x_jacobian(nullptr)
  , quadrature_weights(
      other.matrix_info == nullptr ?
        nullptr :
        mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
  , cell(numbers::invalid_unsigned_int)
  , cell_type(internal::MatrixFreeFunctions::general)
  , is_interior_face(other.is_interior_face)
  , dof_access_index(other.dof_access_index)
  , dof_values_initialized(false)
  , values_quad_initialized(false)
  , gradients_quad_initialized(false)
  , hessians_quad_initialized(false)
  , values_quad_submitted(false)
  , gradients_quad_submitted(false)
  , first_selected_component(other.first_selected_component)
{
  set_data_pointers();
 
  // Create deep copy of mapped geometry for use in parallel...
  if (other.mapped_geometry.get() != nullptr)
    {
      mapped_geometry = std::make_shared<
        internal::MatrixFreeFunctions::
          MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
        other.mapped_geometry->get_fe_values().get_mapping(),
        other.mapped_geometry->get_quadrature(),
        other.mapped_geometry->get_fe_values().get_update_flags());
      mapping_data = &mapped_geometry->get_data_storage();
      cell         = 0;
 
      jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
      J_value  = mapped_geometry->get_data_storage().JxW_values.begin();
    }
}
 
 
 
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)
{
  AssertDimension(quad_no, other.quad_no);
  AssertDimension(n_fe_components, other.n_fe_components);
  AssertDimension(active_fe_index, other.active_fe_index);
  AssertDimension(active_quad_index, other.active_quad_index);
  AssertDimension(first_selected_component, other.first_selected_component);
 
  // release old memory
  if (matrix_info == nullptr)
    {
      delete data;
      delete scratch_data_array;
    }
  else
    {
      matrix_info->release_scratch_data(scratch_data_array);
    }
 
  matrix_info  = other.matrix_info;
  dof_info     = other.dof_info;
  mapping_data = other.mapping_data;
  if (other.matrix_info == nullptr)
    {
      data = new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
        *other.data);
      scratch_data_array = new AlignedVector<VectorizedArrayType>();
    }
  else
    {
      data               = other.data;
      scratch_data_array = matrix_info->acquire_scratch_data();
    }
  set_data_pointers();
 
  quadrature_weights =
    (mapping_data != nullptr ?
       mapping_data->descriptor[active_quad_index].quadrature_weights.begin() :
       nullptr);
  cell             = numbers::invalid_unsigned_int;
  cell_type        = internal::MatrixFreeFunctions::general;
  is_interior_face = other.is_interior_face;
  dof_access_index = other.dof_access_index;
 
  // Create deep copy of mapped geometry for use in parallel...
  if (other.mapped_geometry.get() != nullptr)
    {
      mapped_geometry = std::make_shared<
        internal::MatrixFreeFunctions::
          MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
        other.mapped_geometry->get_fe_values().get_mapping(),
        other.mapped_geometry->get_quadrature(),
        other.mapped_geometry->get_fe_values().get_update_flags());
      cell         = 0;
      mapping_data = &mapped_geometry->get_data_storage();
      jacobian     = mapped_geometry->get_data_storage().jacobians[0].begin();
      J_value      = mapped_geometry->get_data_storage().JxW_values.begin();
    }
 
  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_info != nullptr)
    {
      try
        {
          matrix_info->release_scratch_data(scratch_data_array);
        }
      catch (...)
        {}
    }
  else
    {
      delete scratch_data_array;
      delete data;
      data = nullptr;
    }
  scratch_data_array = nullptr;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  set_data_pointers()
{
  Assert(scratch_data_array != nullptr, ExcInternalError());
 
  const unsigned int tensor_dofs_per_component =
    Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
  const unsigned int dofs_per_component =
    this->data->dofs_per_component_on_cell;
  const unsigned int n_quadrature_points =
    is_face ? this->data->n_q_points_face : this->data->n_q_points;
 
  const unsigned int shift =
    std::max(tensor_dofs_per_component + 1, dofs_per_component) *
      n_components_ * 3 +
    2 * n_quadrature_points;
  const unsigned int allocated_size =
    shift + n_components_ * dofs_per_component +
    (n_components_ * (dim * dim + 2 * dim + 1) * n_quadrature_points);
  scratch_data_array->resize_fast(allocated_size);
 
  // set the pointers to the correct position in the data array
  for (unsigned int c = 0; c < n_components_; ++c)
    {
      this->values_dofs[c] =
        scratch_data_array->begin() + c * dofs_per_component;
      this->values_quad[c] = scratch_data_array->begin() +
                             n_components * dofs_per_component +
                             c * n_quadrature_points;
      for (unsigned int d = 0; d < dim; ++d)
        this->gradients_quad[c][d] =
          scratch_data_array->begin() +
          n_components * (dofs_per_component + n_quadrature_points) +
          (c * dim + d) * n_quadrature_points;
      for (unsigned int d = 0; d < (dim * dim + dim) / 2; ++d)
        this->hessians_quad[c][d] =
          scratch_data_array->begin() +
          n_components *
            ((dim + 1) * n_quadrature_points + dofs_per_component) +
          (c * (dim * dim + dim) + d) * n_quadrature_points;
    }
  scratch_data =
    scratch_data_array->begin() + n_components_ * dofs_per_component +
    (n_components_ * (dim * dim + 2 * dim + 1) * n_quadrature_points);
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline unsigned int
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_mapping_data_index_offset() const
{
  if (matrix_info == nullptr)
    return 0;
  else
    {
      AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
      return this->mapping_data->data_index_offsets[cell];
    }
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline internal::MatrixFreeFunctions::GeometryType
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_cell_type() const
{
  Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  return cell_type;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_shape_info() const
{
  Assert(data != nullptr, ExcInternalError());
  return *data;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  fill_JxW_values(AlignedVector<VectorizedArrayType> &JxW_values) const
{
  AssertDimension(JxW_values.size(), n_quadrature_points);
  Assert(J_value != nullptr, ExcNotInitialized());
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      VectorizedArrayType J = J_value[0];
      for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
        JxW_values[q] = J * this->quadrature_weights[q];
    }
  else
    for (unsigned int q = 0; q < n_quadrature_points; ++q)
      JxW_values[q] = J_value[q];
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
                             FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_normal_vector(const unsigned int q_index) const
{
  AssertIndexRange(q_index, n_quadrature_points);
  Assert(normal_vectors != nullptr, ExcMessage("Did not call reinit()!"));
  if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
    return normal_vectors[0];
  else
    return normal_vectors[q_index];
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
                             FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::JxW(
  const unsigned int q_index) const
{
  AssertIndexRange(q_index, n_quadrature_points);
  Assert(J_value != nullptr, ExcNotInitialized());
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      Assert(this->quadrature_weights != nullptr, ExcInternalError());
      return J_value[0] * this->quadrature_weights[q_index];
    }
  else
    return J_value[q_index];
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline Tensor<2, dim, VectorizedArrayType>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  inverse_jacobian(const unsigned int q_index) const
{
  AssertIndexRange(q_index, n_quadrature_points);
  Assert(this->jacobian != nullptr, ExcNotImplemented());
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    return jacobian[0];
  else
    return jacobian[q_index];
}
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
std::array<unsigned int, VectorizedArrayType::size()>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  get_cell_ids() const
{
  const unsigned int                n_lanes = VectorizedArrayType::size();
  std::array<unsigned int, n_lanes> cells;
 
  // initialize array
  for (unsigned int i = 0; i < n_lanes; ++i)
    cells[i] = numbers::invalid_unsigned_int;
 
  if ((is_face == false) ||
      (is_face &&
       this->dof_access_index ==
         internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
       this->is_interior_face))
    {
      // cell or interior face face (element-centric loop)
      for (unsigned int i = 0; i < n_lanes; ++i)
        cells[i] = cell * n_lanes + i;
    }
  else if (is_face &&
           this->dof_access_index ==
             internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
           this->is_interior_face == false)
    {
      // exterior face (element-centric loop): 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 = this->cell * n_lanes + i;
          // compute face ID
          unsigned int face_index =
            this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
                                                                  this->face_no,
                                                                  i);
 
          if (face_index == numbers::invalid_unsigned_int)
            continue; // invalid face ID: no neighbor on boundary
 
          // get cell ID on both sides of face
          auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes)
                          .cells_interior[face_index % n_lanes];
          auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
                          .cells_exterior[face_index % n_lanes];
 
          // compare the IDs with the given cell ID
          if (cell_m == cell_this)
            cells[i] = cell_p; // neighbor has the other ID
          else if (cell_p == cell_this)
            cells[i] = cell_m;
        }
    }
  else if (is_face)
    {
      // face-centric faces
      const unsigned int *cells_ =
        is_interior_face ?
          &this->matrix_info->get_face_info(cell).cells_interior[0] :
          &this->matrix_info->get_face_info(cell).cells_exterior[0];
      for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
        if (cells_[i] != numbers::invalid_unsigned_int)
          cells[i] = cells_[i];
    }
 
  return cells;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline VectorizedArrayType
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const
{
  Assert(matrix_info != nullptr, ExcNotImplemented());
  AssertDimension(array.size(),
                  matrix_info->get_task_info().cell_partition_data.back());
 
  // 1) collect ids of cell
  const auto cells = this->get_cell_ids();
 
  // 2) actually gather values
  VectorizedArrayType out = make_vectorized_array<Number>(Number(1.));
  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
    if (cells[i] != numbers::invalid_unsigned_int)
      out[i] = array[cells[i] / VectorizedArrayType::size()]
                    [cells[i] % VectorizedArrayType::size()];
  return out;
}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <typename T>
inline std::array<T, VectorizedArrayType::size()>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
                   &array) const
{
  Assert(matrix_info != nullptr, ExcNotImplemented());
  AssertDimension(array.size(),
                  matrix_info->get_task_info().cell_partition_data.back());
 
  // 1) collect ids of cell
  const auto cells = this->get_cell_ids();
 
  // 2) actually gather values
  std::array<T, VectorizedArrayType::size()> out;
  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
    if (cells[i] != numbers::invalid_unsigned_int)
      out[i] = array[cells[i] / VectorizedArrayType::size()]
                    [cells[i] % VectorizedArrayType::size()];
  return out;
}
 
 
 
namespace internal
{
  // 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 VectorType::BlockType;
 
    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<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];
    }
  };
} // 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,
                       VectorType *           src[],
                       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 (matrix_info == nullptr)
    {
      read_write_operation_global(operation, src);
      return;
    }
 
  Assert(dof_info != nullptr, ExcNotInitialized());
  Assert(matrix_info->indices_initialized() == true, ExcNotInitialized());
  if (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], *dof_info);
      }
  else
    {
      internal::check_vector_compatibility(*src[0], *dof_info);
    }
 
  // Case 2: contiguous indices which use reduced storage of indices and can
  // use vectorized load/store operations -> go to separate function
  AssertIndexRange(cell,
                   dof_info->index_storage_variants[dof_access_index].size());
  if (dof_info->index_storage_variants
        [is_face ? dof_access_index :
                   internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
        [cell] >=
      internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
    {
      read_write_operation_contiguous(operation, src, mask);
      return;
    }
 
  // Case 3: standard operation with one index per degree of freedom -> go on
  // here
  constexpr unsigned int n_lanes = VectorizedArrayType::size();
  Assert(mask.count() == n_lanes,
         ExcNotImplemented("Masking currently not implemented for "
                           "non-contiguous DoF storage"));
 
  std::integral_constant<bool,
                         internal::is_vectorizable<VectorType, Number>::value>
    vector_selector;
 
  const unsigned int dofs_per_component =
    this->data->dofs_per_component_on_cell;
  if (dof_info->index_storage_variants
        [is_face ? dof_access_index :
                   internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
        [cell] ==
      internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved)
    {
      const unsigned int *dof_indices =
        dof_info->dof_indices_interleaved.data() +
        dof_info->row_starts[cell * n_fe_components * n_lanes].first +
        dof_info->component_dof_indices_offset[active_fe_index]
                                              [first_selected_component] *
          n_lanes;
      if (n_components == 1 || 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;
    }
 
  const unsigned int *  dof_indices[n_lanes];
  VectorizedArrayType **values_dofs =
    const_cast<VectorizedArrayType **>(&this->values_dofs[0]);
 
  // Assign the appropriate cell ids for face/cell case and get the pointers
  // to the dof indices of the cells on all lanes
  unsigned int        cells_copied[n_lanes];
  const unsigned int *cells;
  unsigned int        n_vectorization_actual =
    dof_info->n_vectorization_lanes_filled[dof_access_index][cell];
  bool has_constraints = false;
  if (is_face)
    {
      if (dof_access_index ==
          internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
        for (unsigned int v = 0; v < n_vectorization_actual; ++v)
          cells_copied[v] = cell * VectorizedArrayType::size() + v;
      cells = dof_access_index ==
                  internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ?
                &cells_copied[0] :
                (is_interior_face ?
                   &this->matrix_info->get_face_info(cell).cells_interior[0] :
                   &this->matrix_info->get_face_info(cell).cells_exterior[0]);
      for (unsigned int v = 0; v < n_vectorization_actual; ++v)
        {
          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] * n_fe_components +
                                  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].second != my_index_start[0].second)
            has_constraints = true;
 
          dof_indices[v] =
            dof_info->dof_indices.data() + my_index_start[0].first;
        }
      for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
        dof_indices[v] = nullptr;
    }
  else
    {
      AssertIndexRange((cell + 1) * n_lanes * n_fe_components,
                       dof_info->row_starts.size());
      const unsigned int n_components_read =
        n_fe_components > 1 ? n_components : 1;
      for (unsigned int v = 0; v < n_vectorization_actual; ++v)
        {
          const std::pair<unsigned int, unsigned int> *my_index_start =
            &dof_info->row_starts[(cell * n_lanes + v) * n_fe_components +
                                  first_selected_component];
          if (my_index_start[n_components_read].second !=
              my_index_start[0].second)
            has_constraints = true;
          Assert(my_index_start[n_components_read].first ==
                     my_index_start[0].first ||
                   my_index_start[0].first < dof_info->dof_indices.size(),
                 ExcIndexRange(0,
                               my_index_start[0].first,
                               dof_info->dof_indices.size()));
          dof_indices[v] =
            dof_info->dof_indices.data() + my_index_start[0].first;
        }
      for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
        dof_indices[v] = nullptr;
    }
 
  // Case where we have no constraints throughout the whole cell: Can go
  // through the list of DoFs directly
  if (!has_constraints)
    {
      if (n_vectorization_actual < 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]);
      if (n_components == 1 || n_fe_components == 1)
        {
          for (unsigned int v = 0; v < n_vectorization_actual; ++v)
            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_vectorization_actual; ++v)
              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
  if (n_vectorization_actual < 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]);
  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
    {
      const unsigned int cell_index = is_face ? cells[v] : cell * n_lanes + v;
      const unsigned int cell_dof_index =
        cell_index * n_fe_components + first_selected_component;
      const unsigned int n_components_read =
        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)
        {
          Assert(dof_info->row_starts_plain_indices[cell_index] !=
                   numbers::invalid_unsigned_int,
                 ExcNotInitialized());
          dof_indices[v] =
            dof_info->plain_dof_indices.data() +
            dof_info->component_dof_indices_offset[active_fe_index]
                                                  [first_selected_component] +
            dof_info->row_starts_plain_indices[cell_index];
          next_index_indicators = index_indicators;
        }
 
      if (n_components == 1 || 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 =
                matrix_info->constraint_pool_begin(indicator.second);
              const Number *end_pool =
                matrix_info->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 =
                    matrix_info->constraint_pool_begin(indicator.second);
                  const Number *end_pool =
                    matrix_info->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,
                              VectorType *           src[]) const
{
  Assert(!local_dof_indices.empty(), ExcNotInitialized());
 
  unsigned int index =
    first_selected_component * data->dofs_per_component_on_cell;
  for (unsigned int comp = 0; comp < n_components; ++comp)
    {
      for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
           ++i, ++index)
        {
          operation.process_empty(values_dofs[comp][i]);
          operation.process_dof_global(
            local_dof_indices[data->lexicographic_numbering[index]],
            *src[0],
            values_dofs[comp][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,
    VectorType *                                    src[],
    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 ? dof_access_index :
              internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
  const unsigned int n_lanes = mask.count();
 
  const std::vector<unsigned int> &dof_indices_cont =
    dof_info->dof_indices_contiguous[ind];
 
  // Simple case: We have contiguous storage, so we can simply copy out the
  // data
  if (dof_info->index_storage_variants[ind][cell] ==
        internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
          interleaved_contiguous &&
      n_lanes == VectorizedArrayType::size())
    {
      const unsigned int dof_index =
        dof_indices_cont[cell * VectorizedArrayType::size()] +
        dof_info->component_dof_indices_offset[active_fe_index]
                                              [first_selected_component] *
          VectorizedArrayType::size();
      if (n_components == 1 || n_fe_components == 1)
        for (unsigned int comp = 0; comp < n_components; ++comp)
          operation.process_dofs_vectorized(data->dofs_per_component_on_cell,
                                            dof_index,
                                            *src[comp],
                                            values_dofs[comp],
                                            vector_selector);
      else
        operation.process_dofs_vectorized(data->dofs_per_component_on_cell *
                                            n_components,
                                          dof_index,
                                          *src[0],
                                          values_dofs[0],
                                          vector_selector);
      return;
    }
 
  // 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];
 
  unsigned int dof_indices[VectorizedArrayType::size()];
  for (unsigned int v = 0; v < n_filled_lanes; ++v)
    dof_indices[v] =
      dof_indices_cont[cell * VectorizedArrayType::size() + v] +
      dof_info->component_dof_indices_offset[active_fe_index]
                                            [first_selected_component] *
        dof_info->dof_indices_interleave_strides
          [ind][cell * VectorizedArrayType::size() + v];
 
  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())
    {
      if (dof_info->index_storage_variants[ind][cell] ==
          internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
            contiguous)
        {
          if (n_components == 1 || n_fe_components == 1)
            for (unsigned int comp = 0; comp < n_components; ++comp)
              operation.process_dofs_vectorized_transpose(
                data->dofs_per_component_on_cell,
                dof_indices,
                *src[comp],
                values_dofs[comp],
                vector_selector);
          else
            operation.process_dofs_vectorized_transpose(
              data->dofs_per_component_on_cell * n_components,
              dof_indices,
              *src[0],
              &values_dofs[0][0],
              vector_selector);
        }
      else if (dof_info->index_storage_variants[ind][cell] ==
               internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                 interleaved_contiguous_strided)
        {
          if (n_components == 1 || n_fe_components == 1)
            for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++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 < data->dofs_per_component_on_cell;
                   ++i)
                {
                  operation.process_dof_gather(
                    dof_indices,
                    *src[0],
                    (comp * data->dofs_per_component_on_cell + i) *
                      VectorizedArrayType::size(),
                    values_dofs[comp][i],
                    vector_selector);
                }
        }
      else
        {
          Assert(dof_info->index_storage_variants[ind][cell] ==
                   internal::MatrixFreeFunctions::DoFInfo::
                     IndexStorageVariants::interleaved_contiguous_mixed_strides,
                 ExcNotImplemented());
          const unsigned int *offsets =
            &dof_info->dof_indices_interleave_strides
               [ind][VectorizedArrayType::size() * cell];
          if (n_components == 1 || n_fe_components == 1)
            for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++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 < data->dofs_per_component_on_cell;
                   ++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 < data->dofs_per_component_on_cell; ++i)
          operation.process_empty(values_dofs[comp][i]);
        if (dof_info->index_storage_variants[ind][cell] ==
            internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
              contiguous)
          {
            if (n_components == 1 || n_fe_components == 1)
              {
                for (unsigned int v = 0; v < n_filled_lanes; ++v)
                  if (mask[v] == true)
                    for (unsigned int i = 0;
                         i < data->dofs_per_component_on_cell;
                         ++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 < data->dofs_per_component_on_cell;
                         ++i)
                      operation.process_dof(
                        dof_indices[v] + i +
                          comp * data->dofs_per_component_on_cell,
                        *src[0],
                        values_dofs[comp][i][v]);
              }
          }
        else
          {
            const unsigned int *offsets =
              &dof_info->dof_indices_interleave_strides
                 [ind][VectorizedArrayType::size() * cell];
            for (unsigned int v = 0; v < n_filled_lanes; ++v)
              AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
            if (n_components == 1 || n_fe_components == 1)
              for (unsigned int v = 0; v < n_filled_lanes; ++v)
                {
                  if (mask[v] == true)
                    for (unsigned int i = 0;
                         i < data->dofs_per_component_on_cell;
                         ++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 < data->dofs_per_component_on_cell;
                         ++i)
                      operation.process_dof(
                        dof_indices[v] +
                          (i + comp * data->dofs_per_component_on_cell) *
                            offsets[v],
                        *src[0],
                        values_dofs[comp][i][v]);
              }
          }
      }
}
 
 
 
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)
{
  // select between block vectors and non-block vectors. Note that the number
  // of components is checked in the internal data
  typename internal::BlockVectorSelector<
    VectorType,
    IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
  for (unsigned int d = 0; d < n_components; ++d)
    src_data[d] =
      internal::BlockVectorSelector<VectorType,
                                    IsBlockVector<VectorType>::value>::
        get_vector_component(const_cast<VectorType &>(src), d + first_index);
 
  internal::VectorReader<Number, VectorizedArrayType> reader;
  read_write_operation(reader,
                       src_data,
                       std::bitset<VectorizedArrayType::size()>().flip(),
                       true);
 
#  ifdef DEBUG
  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)
{
  // select between block vectors and non-block vectors. Note that the number
  // of components is checked in the internal data
  typename internal::BlockVectorSelector<
    VectorType,
    IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
  for (unsigned int d = 0; d < n_components; ++d)
    src_data[d] =
      internal::BlockVectorSelector<VectorType,
                                    IsBlockVector<VectorType>::value>::
        get_vector_component(const_cast<VectorType &>(src), d + first_index);
 
  internal::VectorReader<Number, VectorizedArrayType> reader;
  read_write_operation(reader,
                       src_data,
                       std::bitset<VectorizedArrayType::size()>().flip(),
                       false);
 
#  ifdef DEBUG
  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(dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
 
  // select between block vectors and non-block vectors. Note that the number
  // of components is checked in the internal data
  typename internal::BlockVectorSelector<
    VectorType,
    IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
  for (unsigned int d = 0; d < n_components; ++d)
    dst_data[d] = internal::BlockVectorSelector<
      VectorType,
      IsBlockVector<VectorType>::value>::get_vector_component(dst,
                                                              d + first_index);
 
  internal::VectorDistributorLocalToGlobal<Number, VectorizedArrayType>
    distributor;
  read_write_operation(distributor, dst_data, 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(dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
 
  // select between block vectors and non-block vectors. Note that the number
  // of components is checked in the internal data
  typename internal::BlockVectorSelector<
    VectorType,
    IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
  for (unsigned int d = 0; d < n_components; ++d)
    dst_data[d] = internal::BlockVectorSelector<
      VectorType,
      IsBlockVector<VectorType>::value>::get_vector_component(dst,
                                                              d + first_index);
 
  internal::VectorSetter<Number, VectorizedArrayType> setter;
  read_write_operation(setter, dst_data, mask);
}
 
 
 
/*------------------------------ access to data fields ----------------------*/
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline const std::vector<unsigned int> &
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  get_internal_dof_numbering() const
{
  return data->lexicographic_numbering;
}
 
 
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline ArrayView<VectorizedArrayType>
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  get_scratch_data() const
{
  return ArrayView<VectorizedArrayType>(
    const_cast<VectorizedArrayType *>(scratch_data),
    scratch_data_array->end() - scratch_data);
}
 
 
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline const VectorizedArrayType *
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  begin_dof_values() const
{
  return &values_dofs[0][0];
}
 
 
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline VectorizedArrayType *
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  begin_dof_values()
{
#  ifdef DEBUG
  dof_values_initialized = true;
#  endif
  return &values_dofs[0][0];
}
 
 
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline const VectorizedArrayType *
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  begin_values() const
{
#  ifdef DEBUG
  Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized());
#  endif
  return &values_quad[0][0];
}
 
 
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline VectorizedArrayType *
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  begin_values()
{
#  ifdef DEBUG
  values_quad_initialized = true;
  values_quad_submitted   = true;
#  endif
  return &values_quad[0][0];
}
 
 
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline const VectorizedArrayType *
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  begin_gradients() const
{
#  ifdef DEBUG
  Assert(gradients_quad_initialized || gradients_quad_submitted,
         ExcNotInitialized());
#  endif
  return &gradients_quad[0][0][0];
}
 
 
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline VectorizedArrayType *
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  begin_gradients()
{
#  ifdef DEBUG
  gradients_quad_submitted   = true;
  gradients_quad_initialized = true;
#  endif
  return &gradients_quad[0][0][0];
}
 
 
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline const VectorizedArrayType *
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  begin_hessians() const
{
#  ifdef DEBUG
  Assert(hessians_quad_initialized, ExcNotInitialized());
#  endif
  return &hessians_quad[0][0][0];
}
 
 
 
template <int dim,
          int n_components,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
inline VectorizedArrayType *
FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
  begin_hessians()
{
#  ifdef DEBUG
  hessians_quad_initialized = true;
#  endif
  return &hessians_quad[0][0][0];
}
 
 
 
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);
  Tensor<1, n_components_, VectorizedArrayType> return_value;
  for (unsigned int comp = 0; comp < n_components; comp++)
    return_value[comp] = this->values_dofs[comp][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);
  Tensor<1, n_components_, VectorizedArrayType> return_value;
  for (unsigned int comp = 0; comp < n_components; comp++)
    return_value[comp] = this->values_quad[comp][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(jacobian != nullptr, ExcNotInitialized());
 
  Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_out;
 
  // Cartesian cell
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      for (unsigned int comp = 0; comp < n_components; comp++)
        for (unsigned int d = 0; d < dim; ++d)
          grad_out[comp][d] =
            (this->gradients_quad[comp][d][q_point] * jacobian[0][d][d]);
    }
  // cell with general/affine Jacobian
  else
    {
      const Tensor<2, dim, VectorizedArrayType> &jac =
        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][0][q_point];
            for (unsigned int e = 1; e < dim; ++e)
              grad_out[comp][d] +=
                jac[d][e] * this->gradients_quad[comp][e][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(normal_x_jacobian != nullptr, ExcNotInitialized());
 
  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 - 1][q_point] *
                       (this->normal_x_jacobian[0][dim - 1]);
  else
    {
      const unsigned int 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][0][q_point] *
                           this->normal_x_jacobian[index][0];
          for (unsigned int d = 1; d < dim; ++d)
            grad_out[comp] += this->gradients_quad[comp][d][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_quad[1],
                         const unsigned int               q_point,
                         VectorizedArrayType (&tmp)[1][1])
  {
    tmp[0][0] = jac[0][0] * hessians_quad[0][q_point];
  }
 
  template <typename VectorizedArrayType>
  inline void
  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
                         const VectorizedArrayType *const hessians_quad[3],
                         const unsigned int               q_point,
                         VectorizedArrayType (&tmp)[2][2])
  {
    for (unsigned int d = 0; d < 2; ++d)
      {
        tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
                     jac[d][1] * hessians_quad[2][q_point]);
        tmp[1][d] = (jac[d][0] * hessians_quad[2][q_point] +
                     jac[d][1] * hessians_quad[1][q_point]);
      }
  }
 
  template <typename VectorizedArrayType>
  inline void
  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
                         const VectorizedArrayType *const hessians_quad[6],
                         const unsigned int               q_point,
                         VectorizedArrayType (&tmp)[3][3])
  {
    for (unsigned int d = 0; d < 3; ++d)
      {
        tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
                     jac[d][1] * hessians_quad[3][q_point] +
                     jac[d][2] * hessians_quad[4][q_point]);
        tmp[1][d] = (jac[d][0] * hessians_quad[3][q_point] +
                     jac[d][1] * hessians_quad[1][q_point] +
                     jac[d][2] * hessians_quad[5][q_point]);
        tmp[2][d] = (jac[d][0] * hessians_quad[4][q_point] +
                     jac[d][1] * hessians_quad[5][q_point] +
                     jac[d][2] * hessians_quad[2][q_point]);
      }
  }
} // 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
{
  Assert(!is_face, ExcNotImplemented());
#  ifdef DEBUG
  Assert(this->hessians_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
#  endif
  AssertIndexRange(q_point, this->n_quadrature_points);
 
  Assert(jacobian != nullptr, ExcNotImplemented());
  const Tensor<2, dim, VectorizedArrayType> &jac =
    jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
               0 :
               q_point];
 
  Tensor<2, dim, VectorizedArrayType> hessian_out[n_components];
 
  // 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][d] =
              (this->hessians_quad[comp][d][q_point] * jac[d][d] * jac[d][d]);
            switch (dim)
              {
                case 1:
                  break;
                case 2:
                  hessian_out[comp][0][1] =
                    (this->hessians_quad[comp][2][q_point] * jac[0][0] *
                     jac[1][1]);
                  break;
                case 3:
                  hessian_out[comp][0][1] =
                    (this->hessians_quad[comp][3][q_point] * jac[0][0] *
                     jac[1][1]);
                  hessian_out[comp][0][2] =
                    (this->hessians_quad[comp][4][q_point] * jac[0][0] *
                     jac[2][2]);
                  hessian_out[comp][1][2] =
                    (this->hessians_quad[comp][5][q_point] * jac[1][1] *
                     jac[2][2]);
                  break;
                default:
                  Assert(false, ExcNotImplemented());
              }
            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++)
        {
          // 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],
                                           q_point,
                                           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 Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
        &jac_grad =
          mapping_data->jacobian_gradients
            [1 - this->is_interior_face]
            [this->mapping_data->data_index_offsets[this->cell] + 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],
                                           q_point,
                                           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][e][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][f][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 Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>(
    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(jacobian != nullptr, ExcNotImplemented());
  const Tensor<2, dim, VectorizedArrayType> &jac =
    jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
               0 :
               q_point];
 
  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][d][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],
                                           q_point,
                                           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 Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
        &jac_grad =
          mapping_data->jacobian_gradients
            [0][this->mapping_data->data_index_offsets[this->cell] + 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],
                                           q_point,
                                           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][e][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 Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>
    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
  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
  for (unsigned int comp = 0; comp < n_components; comp++)
    this->values_dofs[comp][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)
{
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr, ExcNotInitialized());
#  ifdef DEBUG
  this->values_quad_submitted = true;
#  endif
 
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        this->values_quad[comp][q_point] = val_in[comp] * JxW;
    }
  else
    {
      const VectorizedArrayType JxW = J_value[q_point];
      for (unsigned int comp = 0; comp < n_components; ++comp)
        this->values_quad[comp][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)
{
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr, ExcNotInitialized());
  Assert(this->jacobian != nullptr, ExcNotInitialized());
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point];
      for (unsigned int comp = 0; comp < n_components; comp++)
        for (unsigned int d = 0; d < dim; ++d)
          this->gradients_quad[comp][d][q_point] =
            (grad_in[comp][d] * jacobian[0][d][d] * JxW);
    }
  else
    {
      const Tensor<2, dim, VectorizedArrayType> &jac =
        this->cell_type > internal::MatrixFreeFunctions::affine ?
          jacobian[q_point] :
          jacobian[0];
      const VectorizedArrayType JxW =
        this->cell_type > internal::MatrixFreeFunctions::affine ?
          J_value[q_point] :
          J_value[0] * 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][d][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, ExcNotInitialized());
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
    for (unsigned int comp = 0; comp < n_components; comp++)
      {
        for (unsigned int d = 0; d < dim - 1; ++d)
          this->gradients_quad[comp][d][q_point] = VectorizedArrayType();
        this->gradients_quad[comp][dim - 1][q_point] =
          grad_in[comp] *
          (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
           this->quadrature_weights[q_point]);
      }
  else
    {
      const unsigned int index =
        this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
      for (unsigned int comp = 0; comp < n_components; comp++)
        {
          VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
          if (this->cell_type <= internal::MatrixFreeFunctions::affine)
            factor = factor * this->quadrature_weights[q_point];
          for (unsigned int d = 0; d < dim; ++d)
            this->gradients_quad[comp][d][q_point] =
              factor * this->normal_x_jacobian[index][d];
        }
    }
}
 
 
 
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
{
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
#  ifdef DEBUG
  Assert(this->values_quad_submitted == true,
         internal::ExcAccessToUninitializedField());
#  endif
  Tensor<1, n_components_, VectorizedArrayType> return_value;
  for (unsigned int comp = 0; comp < n_components; ++comp)
    return_value[comp] = this->values_quad[comp][0];
  const unsigned int n_q_points = this->n_quadrature_points;
  for (unsigned int q = 1; q < n_q_points; ++q)
    for (unsigned int comp = 0; comp < n_components; ++comp)
      return_value[comp] += this->values_quad[comp][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> &data_in,
    const unsigned int                                  dof_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no_in,
    const unsigned int fe_degree,
    const unsigned int n_q_points,
    const bool         is_interior_face)
  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
      data_in,
      dof_no,
      first_selected_component,
      quad_no_in,
      fe_degree,
      n_q_points,
      is_interior_face)
{}
 
 
 
template <int dim,
          int n_components_,
          typename Number,
          bool is_face,
          typename VectorizedArrayType>
template <int n_components_other>
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 FEEvaluationBase<dim,
                                            n_components_other,
                                            Number,
                                            is_face,
                                            VectorizedArrayType> *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> &data_in,
    const unsigned int                                  dof_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no_in,
    const unsigned int fe_degree,
    const unsigned int n_q_points,
    const bool         is_interior_face)
  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
      data_in,
      dof_no,
      first_selected_component,
      quad_no_in,
      fe_degree,
      n_q_points,
      is_interior_face)
{}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
template <int n_components_other>
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 FEEvaluationBase<dim,
                                            n_components_other,
                                            Number,
                                            is_face,
                                            VectorizedArrayType> *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[0][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[0][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, ExcNotInitialized());
 
  Tensor<1, dim, VectorizedArrayType> grad_out;
 
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      for (unsigned int d = 0; d < dim; ++d)
        grad_out[d] =
          (this->gradients_quad[0][d][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[0][0][q_point];
          for (unsigned int e = 1; e < dim; ++e)
            grad_out[d] += jac[d][e] * this->gradients_quad[0][e][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[0][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_index)
{
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_index, this->n_quadrature_points);
  Assert(this->J_value != nullptr, ExcNotInitialized());
#  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_index];
      this->values_quad[0][q_index] = val_in * JxW;
    }
  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
    {
      this->values_quad[0][q_index] = val_in * this->J_value[q_index];
    }
}
 
 
 
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_index)
{
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_index, this->n_quadrature_points);
  Assert(this->J_value != nullptr, ExcNotInitialized());
  Assert(this->jacobian != nullptr, ExcNotInitialized());
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArrayType JxW =
        this->J_value[0] * this->quadrature_weights[q_index];
      for (unsigned int d = 0; d < dim; ++d)
        this->gradients_quad[0][d][q_index] =
          (grad_in[d] * this->jacobian[0][d][d] * JxW);
    }
  // general/affine cell type
  else
    {
      const Tensor<2, dim, VectorizedArrayType> &jac =
        this->cell_type > internal::MatrixFreeFunctions::affine ?
          this->jacobian[q_index] :
          this->jacobian[0];
      const VectorizedArrayType JxW =
        this->cell_type > internal::MatrixFreeFunctions::affine ?
          this->J_value[q_index] :
          this->J_value[0] * this->quadrature_weights[q_index];
      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[0][d][q_index] = new_val * JxW;
        }
    }
}
 
 
 
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> &data_in,
    const unsigned int                                  dof_no,
    const unsigned int first_selected_component,
    const unsigned int quad_no_in,
    const unsigned int fe_degree,
    const unsigned int n_q_points,
    const bool         is_interior_face)
  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
      data_in,
      dof_no,
      first_selected_component,
      quad_no_in,
      fe_degree,
      n_q_points,
      is_interior_face)
{}
 
 
 
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
template <int n_components_other>
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 FEEvaluationBase<dim,
                                            n_components_other,
                                            Number,
                                            is_face,
                                            VectorizedArrayType> *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<2, dim, VectorizedArrayType>
                             FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
  get_gradient(const unsigned int q_point) const
{
  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, ExcNotInitialized());
 
  VectorizedArrayType divergence;
 
  // Cartesian cell
  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      divergence =
        (this->gradients_quad[0][0][q_point] * this->jacobian[0][0][0]);
      for (unsigned int d = 1; d < dim; ++d)
        divergence +=
          (this->gradients_quad[d][d][q_point] * this->jacobian[0][d][d]);
    }
  // cell with general/constant Jacobian
  else
    {
      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[0][0][q_point]);
      for (unsigned int e = 1; e < dim; ++e)
        divergence += (jac[0][e] * this->gradients_quad[0][e][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][e][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 Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
  VectorizedArrayType                       symmetrized[(dim * dim + dim) / 2];
  VectorizedArrayType half = make_vectorized_array<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_gradient(const Tensor<2, dim, VectorizedArrayType> grad_in,
                  const unsigned int                        q_point)
{
  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)
{
  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)
{
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr, ExcNotInitialized());
  Assert(this->jacobian != nullptr, ExcNotInitialized());
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  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][d][q_point] = (fac * this->jacobian[0][d][d]);
          for (unsigned int e = d + 1; e < dim; ++e)
            {
              this->gradients_quad[d][e][q_point] = VectorizedArrayType();
              this->gradients_quad[e][d][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][e][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)
{
  // 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
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr, ExcNotInitialized());
  Assert(this->jacobian != nullptr, ExcNotInitialized());
#  ifdef DEBUG
  this->gradients_quad_submitted = true;
#  endif
 
  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][d][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][d][q_point] =
              (value * this->jacobian[0][d][d]);
            this->gradients_quad[d][e][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][d][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> &data_in,
                     const unsigned int                                dof_no,
                     const unsigned int first_selected_component,
                     const unsigned int quad_no_in,
                     const unsigned int fe_degree,
                     const unsigned int n_q_points,
                     const bool         is_interior_face)
  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
      data_in,
      dof_no,
      first_selected_component,
      quad_no_in,
      fe_degree,
      n_q_points,
      is_interior_face)
{}
 
 
 
template <typename Number, bool is_face, typename VectorizedArrayType>
template <int n_components_other>
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 FEEvaluationBase<1,
                                            n_components_other,
                                            Number,
                                            is_face,
                                            VectorizedArrayType> *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[0][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[0][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[0][0][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_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[0][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)
{
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  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[0][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[0][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)
{
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  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[0][0][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_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 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> &data_in,
               const unsigned int                                  fe_no,
               const unsigned int                                  quad_no,
               const unsigned int first_selected_component)
  : BaseClass(data_in,
              fe_no,
              first_selected_component,
              quad_no,
              fe_degree,
              static_n_q_points)
  , 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 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,
              static_cast<
                FEEvaluationBase<dim, 1, Number, false, VectorizedArrayType> *>(
                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,
              static_cast<
                FEEvaluationBase<dim, 1, Number, false, VectorizedArrayType> *>(
                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>
template <int n_components_other>
inline FEEvaluation<dim,
                    fe_degree,
                    n_q_points_1d,
                    n_components_,
                    Number,
                    VectorizedArrayType>::
  FEEvaluation(const FiniteElement<dim> &                   fe,
               const FEEvaluationBase<dim,
                                      n_components_other,
                                      Number,
                                      false,
                                      VectorizedArrayType> &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;
 
#  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_info->n_components();
                 ++no)
              for (unsigned int nf = 0;
                   nf < this->matrix_info->n_base_elements(no);
                   ++nf)
                if (this->matrix_info
                      ->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_info->get_mapping_info().cell_data.size();
                 ++no)
              if (this->matrix_info->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->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->dof_info != nullptr, ExcNotInitialized());
  Assert(this->mapping_data != nullptr, ExcNotInitialized());
  this->cell = cell_index;
  this->cell_type =
    this->matrix_info->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];
 
#  ifdef DEBUG
  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 <typename DoFHandlerType, bool level_dof_access>
inline void
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::
  reinit(
    const TriaIterator<DoFCellAccessor<DoFHandlerType, level_dof_access>> &cell)
{
  Assert(this->matrix_info == 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().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);
}
 
 
 
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_info == 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);
}
 
 
 
template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number,
          typename VectorizedArrayType>
inline Point<dim, VectorizedArrayType>
FEEvaluation<dim,
             fe_degree,
             n_q_points_1d,
             n_components_,
             Number,
             VectorizedArrayType>::quadrature_point(const unsigned int q) const
{
  if (this->matrix_info == nullptr)
    {
      Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
              update_quadrature_points),
             ExcNotInitialized());
    }
  else
    {
      Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
             ExcNotInitialized());
    }
 
  AssertIndexRange(q, n_q_points);
 
  const Point<dim, VectorizedArrayType> *quadrature_points =
    &this->mapping_data->quadrature_points
       [this->mapping_data->quadrature_point_offsets[this->cell]];
 
  // Cartesian/affine mesh: only first vertex of cell is stored, we must
  // compute it through the Jacobian (which is stored in non-inverted and
  // non-transposed form as index '1' in the jacobian field)
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      Assert(this->jacobian != nullptr, ExcNotInitialized());
      Point<dim, VectorizedArrayType> point = quadrature_points[0];
 
      const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
      if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
        for (unsigned int d = 0; d < dim; ++d)
          point[d] += jac[d][d] *
                      static_cast<Number>(
                        this->mapping_data->descriptor[this->active_quad_index]
                          .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->mapping_data
                                        ->descriptor[this->active_quad_index]
                                        .quadrature.point(q)[e]);
      return point;
    }
  else
    return quadrature_points[q];
}
 
 
 
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[0],
           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 VectorizedArrayType
                                              *        values_array,
                                            const bool evaluate_values,
                                            const bool evaluate_gradients,
                                            const bool evaluate_hessians)
{
  SelectEvaluator<
    dim,
    fe_degree,
    n_q_points_1d,
    n_components,
    VectorizedArrayType>::evaluate(*this->data,
                                   const_cast<VectorizedArrayType *>(
                                     values_array),
                                   this->values_quad[0],
                                   this->gradients_quad[0][0],
                                   this->hessians_quad[0][0],
                                   this->scratch_data,
                                   evaluate_values,
                                   evaluate_gradients,
                                   evaluate_hessians);
 
#  ifdef DEBUG
  if (evaluate_values == true)
    this->values_quad_initialized = true;
  if (evaluate_gradients == true)
    this->gradients_quad_initialized = true;
  if (evaluate_hessians == true)
    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)
{
  // 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() call, without reading the vector entries into a
  // separate data field. This saves some operations.
  if (std::is_same<typename VectorType::value_type, Number>::value &&
      this->dof_info->index_storage_variants
          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
          [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
                            IndexStorageVariants::interleaved_contiguous &&
      reinterpret_cast<std::size_t>(
        input_vector.begin() +
        this->dof_info->dof_indices_contiguous
          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
          [this->cell * VectorizedArrayType::size()]) %
          sizeof(VectorizedArrayType) ==
        0)
    {
      const VectorizedArrayType *vec_values =
        reinterpret_cast<const VectorizedArrayType *>(
          input_vector.begin() +
          this->dof_info->dof_indices_contiguous
            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
            [this->cell * VectorizedArrayType::size()] +
          this->dof_info
              ->component_dof_indices_offset[this->active_fe_index]
                                            [this->first_selected_component] *
            VectorizedArrayType::size());
 
      evaluate(vec_values,
               evaluate_values,
               evaluate_gradients,
               evaluate_hessians);
    }
  else
    {
      this->read_dof_values(input_vector);
      evaluate(this->begin_dof_values(),
               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>::integrate(const bool integrate_values,
                                             const bool integrate_gradients)
{
  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
 
#  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)
{
#  ifdef DEBUG
  if (integrate_values == true)
    Assert(this->values_quad_submitted == true,
           internal::ExcAccessToUninitializedField());
  if (integrate_gradients == true)
    Assert(this->gradients_quad_submitted == true,
           internal::ExcAccessToUninitializedField());
#  endif
  Assert(this->matrix_info != nullptr ||
           this->mapped_geometry->is_initialized(),
         ExcNotInitialized());
 
  SelectEvaluator<dim,
                  fe_degree,
                  n_q_points_1d,
                  n_components,
                  VectorizedArrayType>::integrate(*this->data,
                                                  values_array,
                                                  this->values_quad[0],
                                                  this->gradients_quad[0][0],
                                                  this->scratch_data,
                                                  integrate_values,
                                                  integrate_gradients,
                                                  false);
 
#  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)
{
  // 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 integrate() call, without writing temporary results into a
  // separate data field that will later be added into the vector. This saves
  // some operations.
  if (std::is_same<typename VectorType::value_type, Number>::value &&
      this->dof_info->index_storage_variants
          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
          [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
                            IndexStorageVariants::interleaved_contiguous &&
      reinterpret_cast<std::size_t>(
        destination.begin() +
        this->dof_info->dof_indices_contiguous
          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
          [this->cell * VectorizedArrayType::size()]) %
          sizeof(VectorizedArrayType) ==
        0)
    {
      VectorizedArrayType *vec_values = reinterpret_cast<VectorizedArrayType *>(
        destination.begin() +
        this->dof_info->dof_indices_contiguous
          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
          [this->cell * VectorizedArrayType::size()] +
        this->dof_info
            ->component_dof_indices_offset[this->active_fe_index]
                                          [this->first_selected_component] *
          VectorizedArrayType::size());
      SelectEvaluator<dim,
                      fe_degree,
                      n_q_points_1d,
                      n_components,
                      VectorizedArrayType>::integrate(*this->data,
                                                      vec_values,
                                                      this->values_quad[0],
                                                      this
                                                        ->gradients_quad[0][0],
                                                      this->scratch_data,
                                                      integrate_values,
                                                      integrate_gradients,
                                                      true);
    }
  else
    {
      integrate(integrate_values,
                integrate_gradients,
                this->begin_dof_values());
      this->distribute_local_to_global(destination);
    }
}
 
 
 
/*-------------------------- 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)
  : BaseClass(matrix_free,
              dof_no,
              first_selected_component,
              quad_no,
              fe_degree,
              static_n_q_points,
              is_interior_face)
  , 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_face)
{}
 
 
 
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());
  const internal::MatrixFreeFunctions::FaceToCellTopology<
    VectorizedArrayType::size()> &faces =
    this->matrix_info->get_face_info(face_index);
  if (face_index >=
        this->matrix_info->get_task_info().face_partition_data.back() &&
      face_index <
        this->matrix_info->get_task_info().boundary_partition_data.back())
    Assert(this->is_interior_face,
           ExcMessage("Boundary faces do not have a neighbor"));
 
  this->face_no =
    (this->is_interior_face ? faces.interior_face_no : faces.exterior_face_no);
  this->subface_index = faces.subface_index;
  if (this->is_interior_face == true)
    {
      this->subface_index = GeometryInfo<dim>::max_children_per_cell;
      if (faces.face_orientation > 8)
        this->face_orientation = faces.face_orientation - 8;
      else
        this->face_orientation = 0;
    }
  else
    {
      if (faces.face_orientation < 8)
        this->face_orientation = faces.face_orientation;
      else
        this->face_orientation = 0;
    }
 
  this->cell_type = this->matrix_info->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];
 
#  ifdef DEBUG
  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_info->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_info->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_info != nullptr, ExcNotInitialized());
 
  this->cell_type = this->matrix_info->get_mapping_info().cell_type[cell_index];
  this->cell      = cell_index;
  this->face_orientation = 0;
  this->subface_index    = GeometryInfo<dim>::max_children_per_cell;
  this->face_no          = face_number;
  this->dof_access_index =
    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
 
  const unsigned int offsets =
    this->matrix_info->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_info->get_mapping_info()
                     .face_data_by_cells[this->quad_no]
                     .JxW_values.size());
  this->J_value = &this->matrix_info->get_mapping_info()
                     .face_data_by_cells[this->quad_no]
                     .JxW_values[offsets];
  this->normal_vectors = &this->matrix_info->get_mapping_info()
                            .face_data_by_cells[this->quad_no]
                            .normal_vectors[offsets];
  this->jacobian = &this->matrix_info->get_mapping_info()
                      .face_data_by_cells[this->quad_no]
                      .jacobians[!this->is_interior_face][offsets];
  this->normal_x_jacobian =
    &this->matrix_info->get_mapping_info()
       .face_data_by_cells[this->quad_no]
       .normals_times_jacobians[!this->is_interior_face][offsets];
 
#  ifdef DEBUG
  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[0], 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 VectorizedArrayType
                                                  *        values_array,
                                                const bool evaluate_values,
                                                const bool evaluate_gradients)
{
  if (!(evaluate_values + evaluate_gradients))
    return;
 
  internal::FEFaceEvaluationSelector<
    dim,
    fe_degree,
    n_q_points_1d,
    n_components,
    Number,
    VectorizedArrayType>::evaluate(*this->data,
                                   values_array,
                                   this->begin_values(),
                                   this->begin_gradients(),
                                   this->scratch_data,
                                   evaluate_values,
                                   evaluate_gradients,
                                   this->face_no,
                                   this->subface_index,
                                   this->face_orientation,
                                   this->mapping_data
                                     ->descriptor[this->active_fe_index]
                                     .face_orientations);
 
#  ifdef DEBUG
  if (evaluate_values == true)
    this->values_quad_initialized = true;
  if (evaluate_gradients == true)
    this->gradients_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 bool integrate_values,
                                                 const bool integrate_gradients)
{
  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
 
#  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)
{
  if (!(integrate_values + integrate_gradients))
    return;
 
  internal::FEFaceEvaluationSelector<
    dim,
    fe_degree,
    n_q_points_1d,
    n_components,
    Number,
    VectorizedArrayType>::integrate(*this->data,
                                    values_array,
                                    this->begin_values(),
                                    this->begin_gradients(),
                                    this->scratch_data,
                                    integrate_values,
                                    integrate_gradients,
                                    this->face_no,
                                    this->subface_index,
                                    this->face_orientation,
                                    this->mapping_data
                                      ->descriptor[this->active_fe_index]
                                      .face_orientations);
}
 
 
 
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)
{
  static_assert(internal::has_begin<VectorType>::value &&
                  (std::is_same<decltype(std::declval<VectorType>().begin()),
                                double *>::value ||
                   std::is_same<decltype(std::declval<VectorType>().begin()),
                                float *>::value),
                "This function requires a vector type with begin() function "
                "evaluating to a pointer to basic number (float,double). "
                "Use read_dof_values() followed by evaluate() instead.");
 
  if (!internal::FEFaceEvaluationSelector<dim,
                                          fe_degree,
                                          n_q_points_1d,
                                          n_components,
                                          Number,
                                          VectorizedArrayType>::
        gather_evaluate(input_vector.begin(),
                        *this->data,
                        *this->dof_info,
                        this->begin_values(),
                        this->begin_gradients(),
                        this->scratch_data,
                        evaluate_values,
                        evaluate_gradients,
                        this->active_fe_index,
                        this->first_selected_component,
                        this->cell,
                        this->face_no,
                        this->subface_index,
                        this->dof_access_index,
                        this->face_orientation,
                        this->mapping_data->descriptor[this->active_fe_index]
                          .face_orientations))
    {
      this->read_dof_values(input_vector);
      this->evaluate(evaluate_values, evaluate_gradients);
    }
 
#  ifdef DEBUG
  if (evaluate_values == true)
    this->values_quad_initialized = true;
  if (evaluate_gradients == true)
    this->gradients_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)
{
  static_assert(internal::has_begin<VectorType>::value &&
                  (std::is_same<decltype(std::declval<VectorType>().begin()),
                                double *>::value ||
                   std::is_same<decltype(std::declval<VectorType>().begin()),
                                float *>::value),
                "This function requires a vector type with begin() function "
                "evaluating to a pointer to basic number (float,double). "
                "Use integrate() followed by distribute_local_to_global() "
                "instead.");
 
  if (!internal::FEFaceEvaluationSelector<dim,
                                          fe_degree,
                                          n_q_points_1d,
                                          n_components,
                                          Number,
                                          VectorizedArrayType>::
        integrate_scatter(destination.begin(),
                          *this->data,
                          *this->dof_info,
                          this->begin_dof_values(),
                          this->begin_values(),
                          this->begin_gradients(),
                          this->scratch_data,
                          integrate_values,
                          integrate_gradients,
                          this->active_fe_index,
                          this->first_selected_component,
                          this->cell,
                          this->face_no,
                          this->subface_index,
                          this->dof_access_index,
                          this->face_orientation,
                          this->mapping_data->descriptor[this->active_fe_index]
                            .face_orientations))
    {
      // if we arrive here, writing into the destination vector did not succeed
      // because some of the assumptions in integrate_scatter were not
      // fulfilled (e.g. an element or degree that does not support direct
      // writing), so we must do it here
      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 Point<dim, VectorizedArrayType>
FEFaceEvaluation<dim,
                 fe_degree,
                 n_q_points_1d,
                 n_components_,
                 Number,
                 VectorizedArrayType>::quadrature_point(const unsigned int q)
  const
{
  AssertIndexRange(q, n_q_points);
  if (this->dof_access_index < 2)
    {
      Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
             ExcNotImplemented());
      AssertIndexRange(this->cell,
                       this->mapping_data->quadrature_point_offsets.size());
      return this->mapping_data->quadrature_points
        [this->mapping_data->quadrature_point_offsets[this->cell] + q];
    }
  else
    {
      Assert(this->matrix_info->get_mapping_info()
                 .face_data_by_cells[this->quad_no]
                 .quadrature_point_offsets.empty() == false,
             ExcNotImplemented());
      const unsigned int index =
        this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
      AssertIndexRange(index,
                       this->matrix_info->get_mapping_info()
                         .face_data_by_cells[this->quad_no]
                         .quadrature_point_offsets.size());
      return this->matrix_info->get_mapping_info()
        .face_data_by_cells[this->quad_no]
        .quadrature_points[this->matrix_info->get_mapping_info()
                             .face_data_by_cells[this->quad_no]
                             .quadrature_point_offsets[index] +
                           q];
    }
}
 
 
 
/*------------------------- end FEFaceEvaluation ------------------------- */
 
 
#endif // ifndef DOXYGEN
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif