fe_evaluation.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2018 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>

DEAL_II_NAMESPACE_OPEN



// forward declarations
namespace LinearAlgebra
{
  namespace distributed
  {
    template <typename>
    class Vector;
  }
} // namespace LinearAlgebra
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>
class FEEvaluation;


template <int dim, int n_components_, typename Number, bool is_face = false>
class FEEvaluationBase
{
public:
  using number_type = Number;
  using value_type  = Tensor<1, n_components_, VectorizedArray<Number>>;
  using gradient_type =
    Tensor<1, n_components_, Tensor<1, dim, VectorizedArray<Number>>>;
  static constexpr unsigned int dimension    = dim;
  static constexpr unsigned int n_components = n_components_;

  ~FEEvaluationBase();

  DEAL_II_DEPRECATED
  unsigned int
  get_cell_data_number() const;

  unsigned int
  get_mapping_data_index_offset() const;

  internal::MatrixFreeFunctions::GeometryType
  get_cell_type() const;

  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
  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;

  template <typename VectorType>
  void
  set_dof_values(VectorType &dst, const unsigned int first_index = 0) 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, VectorizedArray<Number>>>
  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:

  VectorizedArray<Number>
  get_divergence(const unsigned int q_point) const;

  SymmetricTensor<2, dim, VectorizedArray<Number>>
  get_symmetric_gradient(const unsigned int q_point) const;

  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArray<Number>>
  get_curl(const unsigned int q_point) const;

  void
  submit_divergence(const VectorizedArray<Number> div_in,
                    const unsigned int            q_point);

  void
  submit_symmetric_gradient(
    const SymmetricTensor<2, dim, VectorizedArray<Number>> grad_in,
    const unsigned int                                     q_point);

  void
  submit_curl(
    const Tensor<1, dim == 2 ? 1 : dim, VectorizedArray<Number>> curl_in,
    const unsigned int                                           q_point);

#endif

  value_type
  integrate_value() const;

  VectorizedArray<Number>
  JxW(const unsigned int q_index) const;

  void
  fill_JxW_values(AlignedVector<VectorizedArray<Number>> &JxW_values) const;

  Tensor<2, dim, VectorizedArray<Number>>
  inverse_jacobian(const unsigned int q_index) const;

  Tensor<1, dim, VectorizedArray<Number>>
  get_normal_vector(const unsigned int q_point) const;

  VectorizedArray<Number>
  read_cell_data(const AlignedVector<VectorizedArray<Number>> &array) const;


  const VectorizedArray<Number> *
  begin_dof_values() const;

  VectorizedArray<Number> *
  begin_dof_values();

  const VectorizedArray<Number> *
  begin_values() const;

  VectorizedArray<Number> *
  begin_values();

  const VectorizedArray<Number> *
  begin_gradients() const;

  VectorizedArray<Number> *
  begin_gradients();

  const VectorizedArray<Number> *
  begin_hessians() const;

  VectorizedArray<Number> *
  begin_hessians();

  const std::vector<unsigned int> &
  get_internal_dof_numbering() const;

  ArrayView<VectorizedArray<Number>>
  get_scratch_data() const;


protected:
  FEEvaluationBase(const MatrixFree<dim, Number> &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> *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 bool             apply_constraints = true) const;

  template <typename VectorType, typename VectorOperation>
  void
  read_write_operation_contiguous(const VectorOperation &operation,
                                  VectorType *           vectors[]) const;

  template <typename VectorType, typename VectorOperation>
  void
  read_write_operation_global(const VectorOperation &operation,
                              VectorType *           vectors[]) const;

  AlignedVector<VectorizedArray<Number>> *scratch_data_array;

  VectorizedArray<Number> *scratch_data;

  VectorizedArray<Number> *values_dofs[n_components];

  VectorizedArray<Number> *values_quad[n_components];

  VectorizedArray<Number> *gradients_quad[n_components][dim];

  VectorizedArray<Number> *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> *matrix_info;

  const internal::MatrixFreeFunctions::DoFInfo *dof_info;

  const internal::MatrixFreeFunctions::
    MappingInfoStorage<(is_face ? dim - 1 : dim), dim, Number> *mapping_data;

  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> *data;

  const Tensor<2, dim, VectorizedArray<Number>> *jacobian;

  const VectorizedArray<Number> *J_value;

  const Tensor<1, dim, VectorizedArray<Number>> *normal_vectors;

  const Tensor<1, dim, VectorizedArray<Number>> *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>>
    mapped_geometry;

  const unsigned int first_selected_component;

  mutable std::vector<types::global_dof_index> local_dof_indices;

private:
  void
  set_data_pointers();

  template <int, int, typename, bool>
  friend class FEEvaluationBase;
  template <int, int, int, int, typename>
  friend class FEEvaluation;
};



template <int dim, int n_components_, typename Number, bool is_face>
class FEEvaluationAccess
  : public FEEvaluationBase<dim, n_components_, Number, is_face>
{
public:
  using number_type = Number;
  using value_type  = Tensor<1, n_components_, VectorizedArray<Number>>;
  using gradient_type =
    Tensor<1, n_components_, Tensor<1, dim, VectorizedArray<Number>>>;
  static constexpr unsigned int dimension    = dim;
  static constexpr unsigned int n_components = n_components_;
  using BaseClass = FEEvaluationBase<dim, n_components_, Number, is_face>;

protected:
  FEEvaluationAccess(const MatrixFree<dim, Number> &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> *other);

  FEEvaluationAccess(const FEEvaluationAccess &other);

  FEEvaluationAccess &
  operator=(const FEEvaluationAccess &other);
};



template <int dim, typename Number, bool is_face>
class FEEvaluationAccess<dim, 1, Number, is_face>
  : public FEEvaluationBase<dim, 1, Number, is_face>
{
public:
  using number_type   = Number;
  using value_type    = VectorizedArray<Number>;
  using gradient_type = Tensor<1, dim, VectorizedArray<Number>>;
  static constexpr unsigned int dimension = dim;
  using BaseClass = FEEvaluationBase<dim, 1, Number, is_face>;

  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, VectorizedArray<Number>> 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, VectorizedArray<Number>>
  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> &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> *other);

  FEEvaluationAccess(const FEEvaluationAccess &other);

  FEEvaluationAccess &
  operator=(const FEEvaluationAccess &other);
};



template <int dim, typename Number, bool is_face>
class FEEvaluationAccess<dim, dim, Number, is_face>
  : public FEEvaluationBase<dim, dim, Number, is_face>
{
public:
  using number_type   = Number;
  using value_type    = Tensor<1, dim, VectorizedArray<Number>>;
  using gradient_type = Tensor<2, dim, VectorizedArray<Number>>;
  static constexpr unsigned int dimension    = dim;
  static constexpr unsigned int n_components = dim;
  using BaseClass = FEEvaluationBase<dim, dim, Number, is_face>;

  gradient_type
  get_gradient(const unsigned int q_point) const;

  VectorizedArray<Number>
  get_divergence(const unsigned int q_point) const;

  SymmetricTensor<2, dim, VectorizedArray<Number>>
  get_symmetric_gradient(const unsigned int q_point) const;

  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArray<Number>>
  get_curl(const unsigned int q_point) const;

  Tensor<3, dim, VectorizedArray<Number>>
  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, VectorizedArray<Number>>> grad_in,
    const unsigned int                                            q_point);

  void
  submit_divergence(const VectorizedArray<Number> div_in,
                    const unsigned int            q_point);

  void
  submit_symmetric_gradient(
    const SymmetricTensor<2, dim, VectorizedArray<Number>> grad_in,
    const unsigned int                                     q_point);

  void
  submit_curl(
    const Tensor<1, dim == 2 ? 1 : dim, VectorizedArray<Number>> curl_in,
    const unsigned int                                           q_point);

protected:
  FEEvaluationAccess(const MatrixFree<dim, Number> &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> *other);

  FEEvaluationAccess(const FEEvaluationAccess &other);

  FEEvaluationAccess &
  operator=(const FEEvaluationAccess &other);
};


template <typename Number, bool is_face>
class FEEvaluationAccess<1, 1, Number, is_face>
  : public FEEvaluationBase<1, 1, Number, is_face>
{
public:
  using number_type   = Number;
  using value_type    = VectorizedArray<Number>;
  using gradient_type = Tensor<1, 1, VectorizedArray<Number>>;
  static constexpr unsigned int dimension = 1;
  using BaseClass = FEEvaluationBase<1, 1, Number, is_face>;

  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, VectorizedArray<Number>>
  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> &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> *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>
class FEEvaluation
  : public FEEvaluationAccess<dim, n_components_, Number, false>
{
public:
  using BaseClass = FEEvaluationAccess<dim, n_components_, Number, false>;

  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> &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> &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 VectorizedArray<Number> *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,
            VectorizedArray<Number> *values_array);

  template <typename VectorType>
  void
  integrate_scatter(const bool  integrate_values,
                    const bool  integrate_gradients,
                    VectorType &output_vector);

  Point<dim, VectorizedArray<Number>>
  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>
class FEFaceEvaluation
  : public FEEvaluationAccess<dim, n_components_, Number, true>
{
public:
  using BaseClass = FEEvaluationAccess<dim, n_components_, Number, true>;

  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> &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);

  ~FEFaceEvaluation();

  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 VectorizedArray<Number> *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,
            VectorizedArray<Number> *values_array);

  template <typename VectorType>
  void
  integrate_scatter(const bool  integrate_values,
                    const bool  integrate_gradients,
                    VectorType &output_vector);

  Point<dim, VectorizedArray<Number>>
  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;

protected:
  void
  adjust_for_face_orientation(const bool integrate,
                              const bool values,
                              const bool gradients);
};



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>
inline FEEvaluationBase<dim, n_components_, Number, is_face>::FEEvaluationBase(
  const MatrixFree<dim, Number> &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>::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_size_info().vectorization_length,
                  VectorizedArray<Number>::n_array_elements);
  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 ||
      (int)n_components_ <=
        (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 " +
      Utilities::to_string(n_components) +
      " components. However, "
      "the current base element has only " +
      Utilities::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 " +
      Utilities::to_string(
        first_selected_component -
        dof_info->start_components
          [dof_info->component_to_base_index[first_selected_component]]) +
      " (global index " + Utilities::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>
template <int n_components_other>
inline FEEvaluationBase<dim, n_components_, Number, is_face>::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> *other)
  : scratch_data_array(new AlignedVector<VectorizedArray<Number>>())
  , 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<VectorizedArray<Number>>(
    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>>(
      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>
inline FEEvaluationBase<dim, n_components_, Number, is_face>::FEEvaluationBase(
  const FEEvaluationBase<dim, n_components_, Number, is_face> &other)
  : scratch_data_array(other.matrix_info == nullptr ?
                         new AlignedVector<VectorizedArray<Number>>() :
                         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<VectorizedArray<Number>>(
          *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.reset(
        new internal::MatrixFreeFunctions::MappingDataOnTheFly<dim, Number>(
          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>
inline FEEvaluationBase<dim, n_components_, Number, is_face> &
FEEvaluationBase<dim, n_components_, Number, is_face>::
operator=(const FEEvaluationBase<dim, n_components_, Number, is_face> &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<VectorizedArray<Number>>(
          *other.data);
      scratch_data_array = new AlignedVector<VectorizedArray<Number>>();
    }
  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.reset(
        new internal::MatrixFreeFunctions::MappingDataOnTheFly<dim, Number>(
          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>
inline FEEvaluationBase<dim, n_components_, Number, is_face>::
  ~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>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face>::set_data_pointers()
{
  Assert(scratch_data_array != nullptr, ExcInternalError());

  const unsigned int tensor_dofs_per_component =
    Utilities::fixed_power<dim>(this->data->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>
inline unsigned int
FEEvaluationBase<dim, n_components_, Number, is_face>::get_cell_data_number()
  const
{
  return get_mapping_data_index_offset();
}



template <int dim, int n_components_, typename Number, bool is_face>
inline unsigned int
FEEvaluationBase<dim, n_components_, Number, is_face>::
  get_mapping_data_index_offset() const
{
  if (matrix_info == 0)
    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>
inline internal::MatrixFreeFunctions::GeometryType
FEEvaluationBase<dim, n_components_, Number, is_face>::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>
inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
FEEvaluationBase<dim, n_components_, Number, is_face>::get_shape_info() const
{
  Assert(data != nullptr, ExcInternalError());
  return *data;
}



template <int dim, int n_components_, typename Number, bool is_face>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face>::fill_JxW_values(
  AlignedVector<VectorizedArray<Number>> &JxW_values) const
{
  AssertDimension(JxW_values.size(), n_quadrature_points);
  Assert(J_value != nullptr, ExcNotInitialized());
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      VectorizedArray<Number> 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>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArray<Number>>
                             FEEvaluationBase<dim, n_components_, Number, is_face>::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>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             FEEvaluationBase<dim, n_components_, Number, is_face>::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>
inline Tensor<2, dim, VectorizedArray<Number>>
FEEvaluationBase<dim, n_components_, Number, is_face>::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>
inline VectorizedArray<Number>
FEEvaluationBase<dim, n_components_, Number, is_face>::read_cell_data(
  const AlignedVector<VectorizedArray<Number>> &array) const
{
  Assert(matrix_info != nullptr, ExcNotImplemented());
  AssertDimension(array.size(),
                  matrix_info->get_task_info().cell_partition_data.back());
  if (is_face)
    {
      VectorizedArray<Number> out = make_vectorized_array<Number>(Number(1.));
      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 < VectorizedArray<Number>::n_array_elements;
           ++i)
        if (cells[i] != numbers::invalid_unsigned_int)
          out[i] = array[cells[i] / VectorizedArray<Number>::n_array_elements]
                        [cells[i] % VectorizedArray<Number>::n_array_elements];
      return out;
    }
  else
    return array[cell];
}



namespace internal
{
  // access to generic vectors that have operator ().
  template <typename VectorType>
  inline typename VectorType::value_type &
  vector_access(VectorType &vec, const unsigned int entry)
  {
    return vec(entry);
  }



  // access to distributed MPI vectors that have a local_element(uint)
  // method to access data in local index space, which is what we use in
  // DoFInfo and hence in read_dof_values etc.
  template <typename Number>
  inline Number &
  vector_access(LinearAlgebra::distributed::Vector<Number> &vec,
                const unsigned int                          entry)
  {
    return vec.local_element(entry);
  }



  // this is to make sure that the parallel partitioning in the
  // LinearAlgebra::distributed::Vector is really the same as stored in
  // MatrixFree
  template <typename VectorType>
  inline void
  check_vector_compatibility(
    const VectorType &                            vec,
    const internal::MatrixFreeFunctions::DoFInfo &dof_info)
  {
    (void)vec;
    (void)dof_info;

    AssertDimension(vec.size(), dof_info.vector_partitioner->size());
  }

  template <typename Number>
  inline void
  check_vector_compatibility(
    const LinearAlgebra::distributed::Vector<Number> &vec,
    const internal::MatrixFreeFunctions::DoFInfo &    dof_info)
  {
    (void)vec;
    (void)dof_info;
    Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
           ExcMessage(
             "The parallel layout of the given vector is not "
             "compatible with the parallel partitioning in MatrixFree. "
             "Use MatrixFree::initialize_dof_vector to get a "
             "compatible vector."));
  }

  // A class to use the same code to read from and write to vector
  template <typename Number>
  struct VectorReader
  {
    template <typename VectorType>
    void
    process_dof(const unsigned int index, VectorType &vec, Number &res) const
    {
      res = vector_access(vec, index);
    }

    template <typename VectorType>
    void
    process_dofs_vectorized_transpose(const unsigned int       dofs_per_cell,
                                      const unsigned int *     dof_indices,
                                      VectorType &             vec,
                                      VectorizedArray<Number> *dof_values,
                                      std::integral_constant<bool, true>) const
    {
      ::vectorized_load_and_transpose(dofs_per_cell,
                                            vec.begin(),
                                            dof_indices,
                                            dof_values);
    }


    template <typename VectorType>
    void
    process_dofs_vectorized_transpose(const unsigned int       dofs_per_cell,
                                      const unsigned int *     dof_indices,
                                      VectorType &             vec,
                                      VectorizedArray<Number> *dof_values,
                                      std::integral_constant<bool, false>) const
    {
      for (unsigned int d = 0; d < dofs_per_cell; ++d)
        for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
             ++v)
          dof_values[d][v] = vector_access(vec, dof_indices[v] + d);
    }

    // variant where VectorType::value_type is the same as Number -> can call
    // gather
    template <typename VectorType>
    void
    process_dof_gather(const unsigned int *     indices,
                       VectorType &             vec,
                       const unsigned int       constant_offset,
                       VectorizedArray<Number> &res,
                       std::integral_constant<bool, true>) const
    {
      res.gather(vec.begin() + constant_offset, indices);
    }

    // variant where VectorType::value_type is not the same as Number -> must
    // manually load the data
    template <typename VectorType>
    void
    process_dof_gather(const unsigned int *     indices,
                       VectorType &             vec,
                       const unsigned int       constant_offset,
                       VectorizedArray<Number> &res,
                       std::integral_constant<bool, false>) const
    {
      for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
           ++v)
        res[v] = vector_access(vec, indices[v] + constant_offset);
    }

    template <typename VectorType>
    void
    process_dof_global(const types::global_dof_index index,
                       VectorType &                  vec,
                       Number &                      res) const
    {
      res = const_cast<const VectorType &>(vec)(index);
    }

    void
    pre_constraints(const Number &, Number &res) const
    {
      res = Number();
    }

    template <typename VectorType>
    void
    process_constraint(const unsigned int index,
                       const Number       weight,
                       VectorType &       vec,
                       Number &           res) const
    {
      res += weight * vector_access(vec, index);
    }

    void
    post_constraints(const Number &sum, Number &write_pos) const
    {
      write_pos = sum;
    }

    void
    process_empty(VectorizedArray<Number> &res) const
    {
      res = VectorizedArray<Number>();
    }
  };

  // A class to use the same code to read from and write to vector
  template <typename Number>
  struct VectorDistributorLocalToGlobal
  {
    template <typename VectorType>
    void
    process_dof(const unsigned int index, VectorType &vec, Number &res) const
    {
      vector_access(vec, index) += res;
    }

    template <typename VectorType>
    void
    process_dofs_vectorized_transpose(const unsigned int       dofs_per_cell,
                                      const unsigned int *     dof_indices,
                                      VectorType &             vec,
                                      VectorizedArray<Number> *dof_values,
                                      std::integral_constant<bool, true>) const
    {
      vectorized_transpose_and_store(
        true, dofs_per_cell, dof_values, dof_indices, vec.begin());
    }

    template <typename VectorType>
    void
    process_dofs_vectorized_transpose(const unsigned int       dofs_per_cell,
                                      const unsigned int *     dof_indices,
                                      VectorType &             vec,
                                      VectorizedArray<Number> *dof_values,
                                      std::integral_constant<bool, false>) const
    {
      for (unsigned int d = 0; d < dofs_per_cell; ++d)
        for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
             ++v)
          vector_access(vec, dof_indices[v] + d) += dof_values[d][v];
    }

    // variant where VectorType::value_type is the same as Number -> can call
    // scatter
    template <typename VectorType>
    void
    process_dof_gather(const unsigned int *     indices,
                       VectorType &             vec,
                       const unsigned int       constant_offset,
                       VectorizedArray<Number> &res,
                       std::integral_constant<bool, true>) const
    {
#  if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3
      for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
           ++v)
        vector_access(vec, indices[v] + constant_offset) += res[v];
#  else
      // only use gather in case there is also scatter.
      VectorizedArray<Number> tmp;
      tmp.gather(vec.begin() + constant_offset, indices);
      tmp += res;
      tmp.scatter(indices, vec.begin() + constant_offset);
#  endif
    }

    // variant where VectorType::value_type is not the same as Number -> must
    // manually append all data
    template <typename VectorType>
    void
    process_dof_gather(const unsigned int *     indices,
                       VectorType &             vec,
                       const unsigned int       constant_offset,
                       VectorizedArray<Number> &res,
                       std::integral_constant<bool, false>) const
    {
      for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
           ++v)
        vector_access(vec, indices[v] + constant_offset) += res[v];
    }

    template <typename VectorType>
    void
    process_dof_global(const types::global_dof_index index,
                       VectorType &                  vec,
                       Number &                      res) const
    {
      vec(index) += res;
    }

    void
    pre_constraints(const Number &input, Number &res) const
    {
      res = input;
    }

    template <typename VectorType>
    void
    process_constraint(const unsigned int index,
                       const Number       weight,
                       VectorType &       vec,
                       Number &           res) const
    {
      vector_access(vec, index) += weight * res;
    }

    void
    post_constraints(const Number &, Number &) const
    {}

    void
    process_empty(VectorizedArray<Number> &) const
    {}
  };


  // A class to use the same code to read from and write to vector
  template <typename Number>
  struct VectorSetter
  {
    template <typename VectorType>
    void
    process_dof(const unsigned int index, VectorType &vec, Number &res) const
    {
      vector_access(vec, index) = res;
    }

    template <typename VectorType>
    void
    process_dofs_vectorized_transpose(const unsigned int       dofs_per_cell,
                                      const unsigned int *     dof_indices,
                                      VectorType &             vec,
                                      VectorizedArray<Number> *dof_values,
                                      std::integral_constant<bool, true>) const
    {
      vectorized_transpose_and_store(
        false, dofs_per_cell, dof_values, dof_indices, vec.begin());
    }

    template <typename VectorType, bool booltype>
    void
    process_dofs_vectorized_transpose(const unsigned int       dofs_per_cell,
                                      const unsigned int *     dof_indices,
                                      VectorType &             vec,
                                      VectorizedArray<Number> *dof_values,
                                      std::integral_constant<bool, false>) const
    {
      for (unsigned int i = 0; i < dofs_per_cell; ++i)
        for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
             ++v)
          vector_access(vec, dof_indices[v] + i) = dof_values[i][v];
    }

    template <typename VectorType>
    void
    process_dof_gather(const unsigned int *     indices,
                       VectorType &             vec,
                       const unsigned int       constant_offset,
                       VectorizedArray<Number> &res,
                       std::integral_constant<bool, true>) const
    {
      res.scatter(indices, vec.begin() + constant_offset);
    }

    template <typename VectorType>
    void
    process_dof_gather(const unsigned int *     indices,
                       VectorType &             vec,
                       const unsigned int       constant_offset,
                       VectorizedArray<Number> &res,
                       std::integral_constant<bool, false>) const
    {
      for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
           ++v)
        vector_access(vec, indices[v] + constant_offset) = res[v];
    }

    template <typename VectorType>
    void
    process_dof_global(const types::global_dof_index index,
                       VectorType &                  vec,
                       Number &                      res) const
    {
      vec(index) = res;
    }

    void
    pre_constraints(const Number &, Number &) const
    {}

    template <typename VectorType>
    void
    process_constraint(const unsigned int,
                       const Number,
                       VectorType &,
                       Number &) const
    {}

    void
    post_constraints(const Number &, Number &) const
    {}

    void
    process_empty(VectorizedArray<Number> &) const
    {}
  };

  // 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>
template <typename VectorType, typename VectorOperation>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face>::read_write_operation(
  const VectorOperation &operation,
  VectorType *           src[],
  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)
      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);
      return;
    }

  // Case 3: standard operation with one index per degree of freedom -> go on
  // here

  constexpr unsigned int n_vectorization =
    VectorizedArray<Number>::n_array_elements;
  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_vectorization].first +
        dof_info->component_dof_indices_offset[active_fe_index]
                                              [first_selected_component] *
          n_vectorization;
      if (n_components == 1 || n_fe_components == 1)
        for (unsigned int i = 0; i < dofs_per_component;
             ++i, dof_indices += n_vectorization)
          for (unsigned int comp = 0; comp < n_components; ++comp)
            operation.process_dof_gather(
              dof_indices,
              *src[comp],
              0,
              values_dofs[comp][i],
              std::integral_constant<
                bool,
                std::is_same<typename VectorType::value_type,
                             Number>::value>());
      else
        for (unsigned int comp = 0; comp < n_components; ++comp)
          for (unsigned int i = 0; i < dofs_per_component;
               ++i, dof_indices += n_vectorization)
            operation.process_dof_gather(
              dof_indices,
              *src[0],
              0,
              values_dofs[comp][i],
              std::integral_constant<
                bool,
                std::is_same<typename VectorType::value_type,
                             Number>::value>());
      return;
    }

  const unsigned int *      dof_indices[n_vectorization];
  VectorizedArray<Number> **values_dofs =
    const_cast<VectorizedArray<Number> **>(&this->values_dofs[0]);

  unsigned int        cells_copied[n_vectorization];
  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 * VectorizedArray<Number>::n_array_elements + 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());
          has_constraints =
            has_constraints &&
            dof_info
                ->row_starts[cells[v] * n_fe_components +
                             first_selected_component + n_components]
                .second != dof_info
                             ->row_starts[cells[v] * n_fe_components +
                                          first_selected_component]
                             .second;
          dof_indices[v] = dof_info->dof_indices.data() +
                           dof_info
                             ->row_starts[cells[v] * n_fe_components +
                                          first_selected_component]
                             .first;
        }
      for (unsigned int v = n_vectorization_actual; v < n_vectorization; ++v)
        dof_indices[v] = nullptr;
    }
  else
    {
      AssertIndexRange((cell + 1) * n_vectorization * 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)
        {
          if (dof_info
                ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                             first_selected_component + n_components_read]
                .second !=
              dof_info
                ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                             first_selected_component]
                .second)
            has_constraints = true;
          Assert(
            dof_info
                  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                               first_selected_component + n_components_read]
                  .first ==
                dof_info
                  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                               first_selected_component]
                  .first ||
              dof_info
                  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                               first_selected_component]
                  .first < dof_info->dof_indices.size(),
            ExcIndexRange(
              0,
              dof_info
                ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                             first_selected_component]
                .first,
              dof_info->dof_indices.size()));
          dof_indices[v] =
            dof_info->dof_indices.data() +
            dof_info
              ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                           first_selected_component]
              .first;
        }
      for (unsigned int v = n_vectorization_actual; v < n_vectorization; ++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_vectorization)
        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_vectorization)
    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)
    {
      unsigned int       index_indicators, next_index_indicators;
      const unsigned int n_components_read =
        n_fe_components > 1 ? n_components : 1;
      if (is_face)
        {
          index_indicators = dof_info
                               ->row_starts[cells[v] * n_fe_components +
                                            first_selected_component]
                               .second;
          next_index_indicators = dof_info
                                    ->row_starts[cells[v] * n_fe_components +
                                                 first_selected_component + 1]
                                    .second;
        }
      else
        {
          index_indicators =
            dof_info
              ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                           first_selected_component]
              .second;
          next_index_indicators =
            dof_info
              ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                           first_selected_component + 1]
              .second;
        }

      if (apply_constraints == false &&
          dof_info
              ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                           first_selected_component]
              .second !=
            dof_info
              ->row_starts[(cell * n_vectorization + v) * n_fe_components +
                           first_selected_component + n_components_read]
              .second)
        {
          Assert(
            dof_info->row_starts_plain_indices[cell * n_vectorization + v] !=
              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] +
            (is_face ?
               dof_info->row_starts_plain_indices[cells[v]] :
               dof_info->row_starts_plain_indices[cell * n_vectorization + v]);
          next_index_indicators = index_indicators;
        }

      if (n_components == 1 || n_fe_components == 1)
        {
          for (unsigned int c = 0; c < n_components; ++c)
            Assert(src[c] != 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."));

          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)
                {
                  if (is_face)
                    next_index_indicators =
                      dof_info
                        ->row_starts[cells[v] * n_fe_components +
                                     first_selected_component + comp + 2]
                        .second;
                  else
                    next_index_indicators =
                      dof_info
                        ->row_starts[(cell * n_vectorization + v) *
                                       n_fe_components +
                                     first_selected_component + comp + 2]
                        .second;
                }
            }
        }
    }
}



template <int dim, int n_components_, typename Number, bool is_face>
template <typename VectorType, typename VectorOperation>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face>::
  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>
template <typename VectorType, typename VectorOperation>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face>::
  read_write_operation_contiguous(const VectorOperation &operation,
                                  VectorType *           src[]) 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,
    std::is_same<typename VectorType::value_type, Number>::value>
                                                               vector_selector;
  const internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex ind =
    is_face ? dof_access_index :
              internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;

  const std::vector<unsigned int> &dof_indices_cont =
    dof_info->dof_indices_contiguous[ind];
  const unsigned int vectorization_populated =
    dof_info->n_vectorization_lanes_filled[ind][this->cell];
  unsigned int dof_indices[VectorizedArray<Number>::n_array_elements];
  for (unsigned int v = 0; v < vectorization_populated; ++v)
    dof_indices[v] =
      dof_indices_cont[cell * VectorizedArray<Number>::n_array_elements + v] +
      dof_info->component_dof_indices_offset[active_fe_index]
                                            [first_selected_component];
  for (unsigned int v = vectorization_populated;
       v < VectorizedArray<Number>::n_array_elements;
       ++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 (vectorization_populated == VectorizedArray<Number>::n_array_elements)
    {
      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
    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 (n_components == 1 || n_fe_components == 1)
          for (unsigned int v = 0; v < vectorization_populated; ++v)
            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 < vectorization_populated; ++v)
            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]);
      }
}



template <int dim, int n_components_, typename Number, bool is_face>
template <typename VectorType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face>::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> reader;
  read_write_operation(reader, src_data, true);

#  ifdef DEBUG
  dof_values_initialized = true;
#  endif
}



template <int dim, int n_components_, typename Number, bool is_face>
template <typename VectorType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face>::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> reader;
  read_write_operation(reader, src_data, false);

#  ifdef DEBUG
  dof_values_initialized = true;
#  endif
}



template <int dim, int n_components_, typename Number, bool is_face>
template <typename VectorType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face>::
  distribute_local_to_global(VectorType &       dst,
                             const unsigned int first_index) const
{
  Assert(dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());

  // 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> distributor;
  read_write_operation(distributor, dst_data);
}



template <int dim, int n_components_, typename Number, bool is_face>
template <typename VectorType>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face>::set_dof_values(
  VectorType &       dst,
  const unsigned int first_index) const
{
  Assert(dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());

  // 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> setter;
  read_write_operation(setter, dst_data);
}



/*------------------------------ access to data fields ----------------------*/

template <int dim, int n_components, typename Number, bool is_face>
inline const std::vector<unsigned int> &
FEEvaluationBase<dim, n_components, Number, is_face>::
  get_internal_dof_numbering() const
{
  return data->lexicographic_numbering;
}



template <int dim, int n_components, typename Number, bool is_face>
inline ArrayView<VectorizedArray<Number>>
FEEvaluationBase<dim, n_components, Number, is_face>::get_scratch_data() const
{
  return ArrayView<VectorizedArray<Number>>(
    const_cast<VectorizedArray<Number> *>(scratch_data),
    scratch_data_array->end() - scratch_data);
}



template <int dim, int n_components, typename Number, bool is_face>
inline const VectorizedArray<Number> *
FEEvaluationBase<dim, n_components, Number, is_face>::begin_dof_values() const
{
  return &values_dofs[0][0];
}



template <int dim, int n_components, typename Number, bool is_face>
inline VectorizedArray<Number> *
FEEvaluationBase<dim, n_components, Number, is_face>::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>
inline const VectorizedArray<Number> *
FEEvaluationBase<dim, n_components, Number, is_face>::begin_values() const
{
  Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized());
  return &values_quad[0][0];
}



template <int dim, int n_components, typename Number, bool is_face>
inline VectorizedArray<Number> *
FEEvaluationBase<dim, n_components, Number, is_face>::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>
inline const VectorizedArray<Number> *
FEEvaluationBase<dim, n_components, Number, is_face>::begin_gradients() const
{
  Assert(gradients_quad_initialized || gradients_quad_submitted,
         ExcNotInitialized());
  return &gradients_quad[0][0][0];
}



template <int dim, int n_components, typename Number, bool is_face>
inline VectorizedArray<Number> *
FEEvaluationBase<dim, n_components, Number, is_face>::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>
inline const VectorizedArray<Number> *
FEEvaluationBase<dim, n_components, Number, is_face>::begin_hessians() const
{
  Assert(hessians_quad_initialized, ExcNotInitialized());
  return &hessians_quad[0][0][0];
}



template <int dim, int n_components, typename Number, bool is_face>
inline VectorizedArray<Number> *
FEEvaluationBase<dim, n_components, Number, is_face>::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>
inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArray<Number>>
                             FEEvaluationBase<dim, n_components_, Number, is_face>::get_dof_value(
  const unsigned int dof) const
{
  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
  Tensor<1, n_components_, VectorizedArray<Number>> 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>
inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArray<Number>>
                             FEEvaluationBase<dim, n_components_, Number, is_face>::get_value(
  const unsigned int q_point) const
{
  Assert(this->values_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);
  Tensor<1, n_components_, VectorizedArray<Number>> 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>
inline DEAL_II_ALWAYS_INLINE
  Tensor<1, n_components_, Tensor<1, dim, VectorizedArray<Number>>>
  FEEvaluationBase<dim, n_components_, Number, is_face>::get_gradient(
    const unsigned int q_point) const
{
  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);

  Assert(jacobian != nullptr, ExcNotInitialized());

  Tensor<1, n_components_, Tensor<1, dim, VectorizedArray<Number>>> 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, VectorizedArray<Number>> &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>
inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArray<Number>>
                             FEEvaluationBase<dim, n_components_, Number, is_face>::get_normal_derivative(
  const unsigned int q_point) const
{
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());

  Assert(normal_x_jacobian != nullptr, ExcNotInitialized());

  Tensor<1, n_components, VectorizedArray<Number>> 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 Number>
  inline void
  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArray<Number>> &jac,
                         const VectorizedArray<Number> *const hessians_quad[1],
                         const unsigned int                   q_point,
                         VectorizedArray<Number> (&tmp)[1][1])
  {
    tmp[0][0] = jac[0][0] * hessians_quad[0][q_point];
  }

  template <typename Number>
  inline void
  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArray<Number>> &jac,
                         const VectorizedArray<Number> *const hessians_quad[3],
                         const unsigned int                   q_point,
                         VectorizedArray<Number> (&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 Number>
  inline void
  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArray<Number>> &jac,
                         const VectorizedArray<Number> *const hessians_quad[6],
                         const unsigned int                   q_point,
                         VectorizedArray<Number> (&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>
inline Tensor<1, n_components_, Tensor<2, dim, VectorizedArray<Number>>>
FEEvaluationBase<dim, n_components_, Number, is_face>::get_hessian(
  const unsigned int q_point) const
{
  Assert(!is_face, ExcNotImplemented());
  Assert(this->hessians_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);

  Assert(jacobian != nullptr, ExcNotImplemented());
  const Tensor<2, dim, VectorizedArray<Number>> &jac =
    jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
               0 :
               q_point];

  Tensor<2, dim, VectorizedArray<Number>> 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
          VectorizedArray<Number> 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, VectorizedArray<Number>>> &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
          VectorizedArray<Number> 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, VectorizedArray<Number>>>(
    hessian_out);
}



template <int dim, int n_components_, typename Number, bool is_face>
inline Tensor<1, n_components_, Tensor<1, dim, VectorizedArray<Number>>>
FEEvaluationBase<dim, n_components_, Number, is_face>::get_hessian_diagonal(
  const unsigned int q_point) const
{
  Assert(!is_face, ExcNotImplemented());
  Assert(this->hessians_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);

  Assert(jacobian != nullptr, ExcNotImplemented());
  const Tensor<2, dim, VectorizedArray<Number>> &jac =
    jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
               0 :
               q_point];

  Tensor<1, n_components_, Tensor<1, dim, VectorizedArray<Number>>> 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
          VectorizedArray<Number> 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, VectorizedArray<Number>>> &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
          VectorizedArray<Number> 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>
inline Tensor<1, n_components_, VectorizedArray<Number>>
FEEvaluationBase<dim, n_components_, Number, is_face>::get_laplacian(
  const unsigned int q_point) const
{
  Assert(is_face == false, ExcNotImplemented());
  Assert(this->hessians_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);

  Tensor<1, n_components_, VectorizedArray<Number>> laplacian_out;
  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArray<Number>>>
    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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face>::submit_dof_value(
  const Tensor<1, n_components_, VectorizedArray<Number>> 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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face>::submit_value(
  const Tensor<1, n_components_, VectorizedArray<Number>> val_in,
  const unsigned int                                      q_point)
{
#  ifdef DEBUG
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->J_value != nullptr, ExcNotInitialized());
  this->values_quad_submitted = true;
#  endif

  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      const VectorizedArray<Number> 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 VectorizedArray<Number> 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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face>::submit_gradient(
  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArray<Number>>>
                     grad_in,
  const unsigned int q_point)
{
#  ifdef DEBUG
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  this->gradients_quad_submitted = true;
  Assert(this->J_value != nullptr, ExcNotInitialized());
  Assert(this->jacobian != nullptr, ExcNotInitialized());
#  endif

  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArray<Number> 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, VectorizedArray<Number>> &jac =
        this->cell_type > internal::MatrixFreeFunctions::affine ?
          jacobian[q_point] :
          jacobian[0];
      const VectorizedArray<Number> 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)
          {
            VectorizedArray<Number> 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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationBase<dim, n_components_, Number, is_face>::submit_normal_derivative(
  const Tensor<1, n_components_, VectorizedArray<Number>> grad_in,
  const unsigned int                                      q_point)
{
#  ifdef DEBUG
  AssertIndexRange(q_point, this->n_quadrature_points);
  this->gradients_quad_submitted = true;
  Assert(this->normal_x_jacobian != nullptr, ExcNotInitialized());
#  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] = VectorizedArray<Number>();
        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++)
        {
          VectorizedArray<Number> 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>
inline Tensor<1, n_components_, VectorizedArray<Number>>
FEEvaluationBase<dim, n_components_, Number, is_face>::integrate_value() const
{
#  ifdef DEBUG
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  Assert(this->values_quad_submitted == true,
         internal::ExcAccessToUninitializedField());
#  endif
  Tensor<1, n_components_, VectorizedArray<Number>> 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>
inline FEEvaluationAccess<dim, n_components_, Number, is_face>::
  FEEvaluationAccess(const MatrixFree<dim, Number> &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>(
      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>
template <int n_components_other>
inline FEEvaluationAccess<dim, n_components_, Number, is_face>::
  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> *other)
  : FEEvaluationBase<dim, n_components_, Number, is_face>(
      mapping,
      fe,
      quadrature,
      update_flags,
      first_selected_component,
      other)
{}



template <int dim, int n_components_, typename Number, bool is_face>
inline FEEvaluationAccess<dim, n_components_, Number, is_face>::
  FEEvaluationAccess(
    const FEEvaluationAccess<dim, n_components_, Number, is_face> &other)
  : FEEvaluationBase<dim, n_components_, Number, is_face>(other)
{}



template <int dim, int n_components_, typename Number, bool is_face>
inline FEEvaluationAccess<dim, n_components_, Number, is_face> &
FEEvaluationAccess<dim, n_components_, Number, is_face>::
operator=(const FEEvaluationAccess<dim, n_components_, Number, is_face> &other)
{
  this->FEEvaluationBase<dim, n_components_, Number, is_face>::operator=(other);
  return *this;
}



/*-------------------- FEEvaluationAccess scalar ----------------------------*/


template <int dim, typename Number, bool is_face>
inline FEEvaluationAccess<dim, 1, Number, is_face>::FEEvaluationAccess(
  const MatrixFree<dim, Number> &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>(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>
template <int n_components_other>
inline FEEvaluationAccess<dim, 1, Number, is_face>::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> *other)
  : FEEvaluationBase<dim, 1, Number, is_face>(mapping,
                                              fe,
                                              quadrature,
                                              update_flags,
                                              first_selected_component,
                                              other)
{}



template <int dim, typename Number, bool is_face>
inline FEEvaluationAccess<dim, 1, Number, is_face>::FEEvaluationAccess(
  const FEEvaluationAccess<dim, 1, Number, is_face> &other)
  : FEEvaluationBase<dim, 1, Number, is_face>(other)
{}



template <int dim, typename Number, bool is_face>
inline FEEvaluationAccess<dim, 1, Number, is_face> &
FEEvaluationAccess<dim, 1, Number, is_face>::
operator=(const FEEvaluationAccess<dim, 1, Number, is_face> &other)
{
  this->FEEvaluationBase<dim, 1, Number, is_face>::operator=(other);
  return *this;
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             FEEvaluationAccess<dim, 1, Number, is_face>::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>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             FEEvaluationAccess<dim, 1, Number, is_face>::get_value(
  const unsigned int q_point) const
{
  Assert(this->values_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);
  return this->values_quad[0][q_point];
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             FEEvaluationAccess<dim, 1, Number, is_face>::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>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArray<Number>>
                             FEEvaluationAccess<dim, 1, Number, is_face>::get_gradient(
  const unsigned int q_point) const
{
  // could use the base class gradient, but that involves too many expensive
  // initialization operations on tensors

  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);

  Assert(this->jacobian != nullptr, ExcNotInitialized());

  Tensor<1, dim, VectorizedArray<Number>> 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, VectorizedArray<Number>> &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>
inline Tensor<2, dim, VectorizedArray<Number>>
FEEvaluationAccess<dim, 1, Number, is_face>::get_hessian(
  const unsigned int q_point) const
{
  return BaseClass::get_hessian(q_point)[0];
}



template <int dim, typename Number, bool is_face>
inline Tensor<1, dim, VectorizedArray<Number>>
FEEvaluationAccess<dim, 1, Number, is_face>::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>
inline VectorizedArray<Number>
FEEvaluationAccess<dim, 1, Number, is_face>::get_laplacian(
  const unsigned int q_point) const
{
  return BaseClass::get_laplacian(q_point)[0];
}



template <int dim, typename Number, bool is_face>
inline void DEAL_II_ALWAYS_INLINE
            FEEvaluationAccess<dim, 1, Number, is_face>::submit_dof_value(
  const VectorizedArray<Number> 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>
inline void DEAL_II_ALWAYS_INLINE
            FEEvaluationAccess<dim, 1, Number, is_face>::submit_value(
  const VectorizedArray<Number> val_in,
  const unsigned int            q_index)
{
#  ifdef DEBUG
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_index, this->n_quadrature_points);
  Assert(this->J_value != nullptr, ExcNotInitialized());
  this->values_quad_submitted = true;
#  endif
  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      const VectorizedArray<Number> 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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, 1, Number, is_face>::submit_value(
  const Tensor<1, 1, VectorizedArray<Number>> val_in,
  const unsigned int                          q_point)
{
  submit_value(val_in[0], q_point);
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, 1, Number, is_face>::submit_normal_derivative(
  const VectorizedArray<Number> grad_in,
  const unsigned int            q_point)
{
  Tensor<1, 1, VectorizedArray<Number>> grad;
  grad[0] = grad_in;
  BaseClass::submit_normal_derivative(grad, q_point);
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, 1, Number, is_face>::submit_gradient(
  const Tensor<1, dim, VectorizedArray<Number>> grad_in,
  const unsigned int                            q_index)
{
#  ifdef DEBUG
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_index, this->n_quadrature_points);
  this->gradients_quad_submitted = true;
  Assert(this->J_value != nullptr, ExcNotInitialized());
  Assert(this->jacobian != nullptr, ExcNotInitialized());
#  endif

  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArray<Number> 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, VectorizedArray<Number>> &jac =
        this->cell_type > internal::MatrixFreeFunctions::affine ?
          this->jacobian[q_index] :
          this->jacobian[0];
      const VectorizedArray<Number> 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)
        {
          VectorizedArray<Number> 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>
inline VectorizedArray<Number>
FEEvaluationAccess<dim, 1, Number, is_face>::integrate_value() const
{
  return BaseClass::integrate_value()[0];
}



/*----------------- FEEvaluationAccess vector-valued ------------------------*/


template <int dim, typename Number, bool is_face>
inline FEEvaluationAccess<dim, dim, Number, is_face>::FEEvaluationAccess(
  const MatrixFree<dim, Number> &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>(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>
template <int n_components_other>
inline FEEvaluationAccess<dim, dim, Number, is_face>::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> *other)
  : FEEvaluationBase<dim, dim, Number, is_face>(mapping,
                                                fe,
                                                quadrature,
                                                update_flags,
                                                first_selected_component,
                                                other)
{}



template <int dim, typename Number, bool is_face>
inline FEEvaluationAccess<dim, dim, Number, is_face>::FEEvaluationAccess(
  const FEEvaluationAccess<dim, dim, Number, is_face> &other)
  : FEEvaluationBase<dim, dim, Number, is_face>(other)
{}



template <int dim, typename Number, bool is_face>
inline FEEvaluationAccess<dim, dim, Number, is_face> &
FEEvaluationAccess<dim, dim, Number, is_face>::
operator=(const FEEvaluationAccess<dim, dim, Number, is_face> &other)
{
  this->FEEvaluationBase<dim, dim, Number, is_face>::operator=(other);
  return *this;
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArray<Number>>
                             FEEvaluationAccess<dim, dim, Number, is_face>::get_gradient(
  const unsigned int q_point) const
{
  return BaseClass::get_gradient(q_point);
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             FEEvaluationAccess<dim, dim, Number, is_face>::get_divergence(
  const unsigned int q_point) const
{
  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);
  Assert(this->jacobian != nullptr, ExcNotInitialized());

  VectorizedArray<Number> 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, VectorizedArray<Number>> &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>
inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, VectorizedArray<Number>>
                             FEEvaluationAccess<dim, dim, Number, is_face>::get_symmetric_gradient(
  const unsigned int q_point) const
{
  // copy from generic function into dim-specialization function
  const Tensor<2, dim, VectorizedArray<Number>> grad = get_gradient(q_point);
  VectorizedArray<Number> symmetrized[(dim * dim + dim) / 2];
  VectorizedArray<Number> 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, VectorizedArray<Number>>(symmetrized);
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE
  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArray<Number>>
  FEEvaluationAccess<dim, dim, Number, is_face>::get_curl(
    const unsigned int q_point) const
{
  // copy from generic function into dim-specialization function
  const Tensor<2, dim, VectorizedArray<Number>> grad = get_gradient(q_point);
  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArray<Number>> 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>
inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArray<Number>>
                             FEEvaluationAccess<dim, dim, Number, is_face>::get_hessian_diagonal(
  const unsigned int q_point) const
{
  return BaseClass::get_hessian_diagonal(q_point);
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Tensor<3, dim, VectorizedArray<Number>>
                             FEEvaluationAccess<dim, dim, Number, is_face>::get_hessian(
  const unsigned int q_point) const
{
  Assert(this->hessians_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);
  return BaseClass::get_hessian(q_point);
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face>::submit_gradient(
  const Tensor<2, dim, VectorizedArray<Number>> grad_in,
  const unsigned int                            q_point)
{
  BaseClass::submit_gradient(grad_in, q_point);
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face>::submit_gradient(
  const Tensor<1, dim, Tensor<1, dim, VectorizedArray<Number>>> grad_in,
  const unsigned int                                            q_point)
{
  BaseClass::submit_gradient(grad_in, q_point);
}



template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face>::submit_divergence(
  const VectorizedArray<Number> div_in,
  const unsigned int            q_point)
{
#  ifdef DEBUG
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  this->gradients_quad_submitted = true;
  Assert(this->J_value != nullptr, ExcNotInitialized());
  Assert(this->jacobian != nullptr, ExcNotInitialized());
#  endif

  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArray<Number> 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] = VectorizedArray<Number>();
              this->gradients_quad[e][d][q_point] = VectorizedArray<Number>();
            }
        }
    }
  else
    {
      const Tensor<2, dim, VectorizedArray<Number>> &jac =
        this->cell_type == internal::MatrixFreeFunctions::general ?
          this->jacobian[q_point] :
          this->jacobian[0];
      const VectorizedArray<Number> 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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face>::submit_symmetric_gradient(
  const SymmetricTensor<2, dim, VectorizedArray<Number>> 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
#  ifdef DEBUG
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  this->gradients_quad_submitted = true;
  Assert(this->J_value != nullptr, ExcNotInitialized());
  Assert(this->jacobian != nullptr, ExcNotInitialized());
#  endif

  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      const VectorizedArray<Number> 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 VectorizedArray<Number> 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 VectorizedArray<Number> 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, VectorizedArray<Number>> &jac =
        this->cell_type == internal::MatrixFreeFunctions::general ?
          this->jacobian[q_point] :
          this->jacobian[0];
      VectorizedArray<Number> 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 VectorizedArray<Number> 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)
          {
            VectorizedArray<Number> 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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<dim, dim, Number, is_face>::submit_curl(
  const Tensor<1, dim == 2 ? 1 : dim, VectorizedArray<Number>> curl,
  const unsigned int                                           q_point)
{
  Tensor<2, dim, VectorizedArray<Number>> 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>
inline FEEvaluationAccess<1, 1, Number, is_face>::FEEvaluationAccess(
  const MatrixFree<1, Number> &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>(data_in,
                                            dof_no,
                                            first_selected_component,
                                            quad_no_in,
                                            fe_degree,
                                            n_q_points,
                                            is_interior_face)
{}



template <typename Number, bool is_face>
template <int n_components_other>
inline FEEvaluationAccess<1, 1, Number, is_face>::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> *other)
  : FEEvaluationBase<1, 1, Number, is_face>(mapping,
                                            fe,
                                            quadrature,
                                            update_flags,
                                            first_selected_component,
                                            other)
{}



template <typename Number, bool is_face>
inline FEEvaluationAccess<1, 1, Number, is_face>::FEEvaluationAccess(
  const FEEvaluationAccess<1, 1, Number, is_face> &other)
  : FEEvaluationBase<1, 1, Number, is_face>(other)
{}



template <typename Number, bool is_face>
inline FEEvaluationAccess<1, 1, Number, is_face> &
FEEvaluationAccess<1, 1, Number, is_face>::
operator=(const FEEvaluationAccess<1, 1, Number, is_face> &other)
{
  this->FEEvaluationBase<1, 1, Number, is_face>::operator=(other);
  return *this;
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             FEEvaluationAccess<1, 1, Number, is_face>::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>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             FEEvaluationAccess<1, 1, Number, is_face>::get_value(
  const unsigned int q_point) const
{
  Assert(this->values_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);
  return this->values_quad[0][q_point];
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Tensor<1, 1, VectorizedArray<Number>>
                             FEEvaluationAccess<1, 1, Number, is_face>::get_gradient(
  const unsigned int q_point) const
{
  // could use the base class gradient, but that involves too many inefficient
  // initialization operations on tensors

  Assert(this->gradients_quad_initialized == true,
         internal::ExcAccessToUninitializedField());
  AssertIndexRange(q_point, this->n_quadrature_points);

  const Tensor<2, 1, VectorizedArray<Number>> &jac =
    this->cell_type == internal::MatrixFreeFunctions::general ?
      this->jacobian[q_point] :
      this->jacobian[0];

  Tensor<1, 1, VectorizedArray<Number>> grad_out;
  grad_out[0] = jac[0][0] * this->gradients_quad[0][0][q_point];

  return grad_out;
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             FEEvaluationAccess<1, 1, Number, is_face>::get_normal_derivative(
  const unsigned int q_point) const
{
  return BaseClass::get_normal_derivative(q_point)[0];
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Tensor<2, 1, VectorizedArray<Number>>
                             FEEvaluationAccess<1, 1, Number, is_face>::get_hessian(
  const unsigned int q_point) const
{
  return BaseClass::get_hessian(q_point)[0];
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Tensor<1, 1, VectorizedArray<Number>>
                             FEEvaluationAccess<1, 1, Number, is_face>::get_hessian_diagonal(
  const unsigned int q_point) const
{
  return BaseClass::get_hessian_diagonal(q_point)[0];
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             FEEvaluationAccess<1, 1, Number, is_face>::get_laplacian(
  const unsigned int q_point) const
{
  return BaseClass::get_laplacian(q_point)[0];
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE void DEAL_II_ALWAYS_INLINE
                                  FEEvaluationAccess<1, 1, Number, is_face>::submit_dof_value(
  const VectorizedArray<Number> 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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face>::submit_value(
  const VectorizedArray<Number> val_in,
  const unsigned int            q_point)
{
#  ifdef DEBUG
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  this->values_quad_submitted = true;
#  endif
  if (this->cell_type == internal::MatrixFreeFunctions::general)
    {
      const VectorizedArray<Number> JxW = this->J_value[q_point];
      this->values_quad[0][q_point]     = val_in * JxW;
    }
  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
    {
      const VectorizedArray<Number> 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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face>::submit_value(
  const Tensor<1, 1, VectorizedArray<Number>> val_in,
  const unsigned int                          q_point)
{
  submit_value(val_in[0], q_point);
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face>::submit_gradient(
  const Tensor<1, 1, VectorizedArray<Number>> grad_in,
  const unsigned int                          q_point)
{
  submit_gradient(grad_in[0], q_point);
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face>::submit_gradient(
  const VectorizedArray<Number> grad_in,
  const unsigned int            q_point)
{
#  ifdef DEBUG
  Assert(this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
  AssertIndexRange(q_point, this->n_quadrature_points);
  this->gradients_quad_submitted = true;
#  endif

  const Tensor<2, 1, VectorizedArray<Number>> &jac =
    this->cell_type == internal::MatrixFreeFunctions::general ?
      this->jacobian[q_point] :
      this->jacobian[0];
  const VectorizedArray<Number> 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>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face>::submit_normal_derivative(
  const VectorizedArray<Number> grad_in,
  const unsigned int            q_point)
{
  Tensor<1, 1, VectorizedArray<Number>> grad;
  grad[0] = grad_in;
  BaseClass::submit_normal_derivative(grad, q_point);
}



template <typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE void
FEEvaluationAccess<1, 1, Number, is_face>::submit_normal_derivative(
  const Tensor<1, 1, VectorizedArray<Number>> grad_in,
  const unsigned int                          q_point)
{
  BaseClass::submit_normal_derivative(grad_in, q_point);
}



template <typename Number, bool is_face>
inline VectorizedArray<Number>
FEEvaluationAccess<1, 1, Number, is_face>::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>
inline FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  FEEvaluation(const MatrixFree<dim, Number> &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>
inline FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  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> *>(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>
inline FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  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> *>(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>
template <int n_components_other>
inline FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  FEEvaluation(const FiniteElement<dim> &                               fe,
               const FEEvaluationBase<dim, n_components_other, Number> &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>
inline FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  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>
inline FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number> &
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
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>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  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->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->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)
                      .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->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->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->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>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::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>
template <typename DoFHandlerType, bool level_dof_access>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::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>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::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>
inline Point<dim, VectorizedArray<Number>>
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  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 unsigned int n_q_points_1d_actual =
    fe_degree == -1 ? this->data->n_q_points_1d : n_q_points_1d;

  // Cartesian mesh: not all quadrature points are stored, only the
  // diagonal. Hence, need to find the tensor product index and retrieve the
  // value from that
  const Point<dim, VectorizedArray<Number>> *quadrature_points =
    &this->mapping_data->quadrature_points
       [this->mapping_data->quadrature_point_offsets[this->cell]];
  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
    {
      Point<dim, VectorizedArray<Number>> point;
      switch (dim)
        {
          case 1:
            return quadrature_points[q];
          case 2:
            point[0] = quadrature_points[q % n_q_points_1d_actual][0];
            point[1] = quadrature_points[q / n_q_points_1d_actual][1];
            return point;
          case 3:
            point[0] = quadrature_points[q % n_q_points_1d_actual][0];
            point[1] = quadrature_points[(q / n_q_points_1d_actual) %
                                         n_q_points_1d_actual][1];
            point[2] = quadrature_points[q / (n_q_points_1d_actual *
                                              n_q_points_1d_actual)][2];
            return point;
          default:
            Assert(false, ExcNotImplemented());
            return point;
        }
    }
  // all other cases: just return the respective data as it is fully stored
  else
    return quadrature_points[q];
}



template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::evaluate(
  const bool evaluate_values,
  const bool evaluate_gradients,
  const bool evaluate_hessians)
{
  Assert(this->dof_values_initialized == true,
         internal::ExcAccessToUninitializedField());
  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>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::evaluate(
  const VectorizedArray<Number> *values_array,
  const bool                     evaluate_values,
  const bool                     evaluate_gradients,
  const bool                     evaluate_hessians)
{
  SelectEvaluator<
    dim,
    fe_degree,
    n_q_points_1d,
    n_components,
    VectorizedArray<Number>>::evaluate(*this->data,
                                       const_cast<VectorizedArray<Number> *>(
                                         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>
template <typename VectorType>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  gather_evaluate(const VectorType &input_vector,
                  const bool        evaluate_values,
                  const bool        evaluate_gradients,
                  const bool        evaluate_hessians)
{
  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>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::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>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::integrate(
  const bool               integrate_values,
  const bool               integrate_gradients,
  VectorizedArray<Number> *values_array)
{
  if (integrate_values == true)
    Assert(this->values_quad_submitted == true,
           internal::ExcAccessToUninitializedField());
  if (integrate_gradients == true)
    Assert(this->gradients_quad_submitted == true,
           internal::ExcAccessToUninitializedField());
  Assert(this->matrix_info != nullptr ||
           this->mapped_geometry->is_initialized(),
         ExcNotInitialized());

  SelectEvaluator<dim,
                  fe_degree,
                  n_q_points_1d,
                  n_components,
                  VectorizedArray<Number>>::integrate(*this->data,
                                                      values_array,
                                                      this->values_quad[0],
                                                      this
                                                        ->gradients_quad[0][0],
                                                      this->scratch_data,
                                                      integrate_values,
                                                      integrate_gradients);

#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif
}



template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number>
template <typename VectorType>
inline void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  integrate_scatter(const bool  integrate_values,
                    const bool  integrate_gradients,
                    VectorType &destination)
{
  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>
inline FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  FEFaceEvaluation(const MatrixFree<dim, Number> &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>
inline FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  ~FEFaceEvaluation()
{}



template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number>
inline void
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::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 unsigned int n_vectors = VectorizedArray<Number>::n_array_elements;
  const internal::MatrixFreeFunctions::FaceToCellTopology<n_vectors> &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->values_quad_submitted = false;

  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>
inline void
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::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"));
  Assert(this->is_interior_face == true,
         ExcMessage(
           "Cell-based FEFaceEvaluation::reinit only possible for the "
           "interior face with second argument to constructor as true"));
  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[0][offsets];
  this->normal_x_jacobian = &this->matrix_info->get_mapping_info()
                               .face_data_by_cells[this->quad_no]
                               .normals_times_jacobians[0][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>
inline void
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>::evaluate(
  const bool evaluate_values,
  const bool evaluate_gradients)
{
  Assert(this->dof_values_initialized, ExcNotInitialized());

  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>
inline void
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>::evaluate(
  const VectorizedArray<Number> *values_array,
  const bool                     evaluate_values,
  const bool                     evaluate_gradients)
{
  if (!(evaluate_values + evaluate_gradients))
    return;

  constexpr unsigned int static_dofs_per_face =
    fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
                     numbers::invalid_unsigned_int;
  const unsigned int dofs_per_face =
    fe_degree > -1 ? static_dofs_per_face :
                     Utilities::pow(this->data->fe_degree + 1, dim - 1);

  // we allocate small amounts of data on the stack to signal the compiler
  // that this temporary data is only needed for the calculations but the
  // final results can be discarded and need not be written back to
  // memory. For large sizes or when the dofs per face is not a compile-time
  // constant, however, we want to go to the heap in the `scratch_data`
  // variable to not risk a stack overflow.
  constexpr unsigned int stack_array_size_threshold = 100;

  VectorizedArray<Number>
                           temp_data[static_dofs_per_face < stack_array_size_threshold ?
                n_components * 2 * static_dofs_per_face :
                1];
  VectorizedArray<Number> *temp1;
  if (static_dofs_per_face < stack_array_size_threshold)
    temp1 = &temp_data[0];
  else
    temp1 = this->scratch_data;

  internal::FEFaceNormalEvaluationImpl<dim,
                                       fe_degree,
                                       n_components,
                                       VectorizedArray<Number>>::
    template interpolate<true, false>(
      *this->data, values_array, temp1, evaluate_gradients, this->face_no);

  const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0;
  if (fe_degree > -1 &&
      this->subface_index >= GeometryInfo<dim>::max_children_per_cell &&
      this->data->element_type <=
        internal::MatrixFreeFunctions::tensor_symmetric)
    internal::FEFaceEvaluationImpl<
      true,
      dim,
      fe_degree,
      n_q_points_1d_actual,
      n_components,
      VectorizedArray<Number>>::evaluate_in_face(*this->data,
                                                 temp1,
                                                 this->begin_values(),
                                                 this->begin_gradients(),
                                                 this->scratch_data +
                                                   2 * n_components *
                                                     dofs_per_face,
                                                 evaluate_values,
                                                 evaluate_gradients,
                                                 this->subface_index);
  else
    internal::FEFaceEvaluationImpl<
      false,
      dim,
      fe_degree,
      n_q_points_1d_actual,
      n_components,
      VectorizedArray<Number>>::evaluate_in_face(*this->data,
                                                 temp1,
                                                 this->begin_values(),
                                                 this->begin_gradients(),
                                                 this->scratch_data +
                                                   2 * n_components *
                                                     dofs_per_face,
                                                 evaluate_values,
                                                 evaluate_gradients,
                                                 this->subface_index);

  if (this->face_orientation)
    adjust_for_face_orientation(false, 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>
inline void
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>::
  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>
inline void
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>::
  integrate(const bool               integrate_values,
            const bool               integrate_gradients,
            VectorizedArray<Number> *values_array)
{
  if (!(integrate_values + integrate_gradients))
    return;

  if (this->face_orientation)
    adjust_for_face_orientation(true, integrate_values, integrate_gradients);

  constexpr unsigned int static_dofs_per_face =
    fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
                     numbers::invalid_unsigned_int;
  const unsigned int dofs_per_face =
    fe_degree > -1 ? static_dofs_per_face :
                     Utilities::pow(this->data->fe_degree + 1, dim - 1);

  constexpr unsigned int stack_array_size_threshold = 100;

  VectorizedArray<Number>
                           temp_data[static_dofs_per_face < stack_array_size_threshold ?
                n_components * 2 * static_dofs_per_face :
                1];
  VectorizedArray<Number> *temp1;
  if (static_dofs_per_face < stack_array_size_threshold)
    temp1 = &temp_data[0];
  else
    temp1 = this->scratch_data;

  const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0;
  if (fe_degree > -1 &&
      this->subface_index >= GeometryInfo<dim - 1>::max_children_per_cell &&
      this->data->element_type <=
        internal::MatrixFreeFunctions::tensor_symmetric)
    internal::FEFaceEvaluationImpl<
      true,
      dim,
      fe_degree,
      n_q_points_1d_actual,
      n_components,
      VectorizedArray<Number>>::integrate_in_face(*this->data,
                                                  temp1,
                                                  this->begin_values(),
                                                  this->begin_gradients(),
                                                  this->scratch_data +
                                                    2 * n_components *
                                                      dofs_per_face,
                                                  integrate_values,
                                                  integrate_gradients,
                                                  this->subface_index);
  else
    internal::FEFaceEvaluationImpl<
      false,
      dim,
      fe_degree,
      n_q_points_1d_actual,
      n_components,
      VectorizedArray<Number>>::integrate_in_face(*this->data,
                                                  temp1,
                                                  this->begin_values(),
                                                  this->begin_gradients(),
                                                  this->scratch_data +
                                                    2 * n_components *
                                                      dofs_per_face,
                                                  integrate_values,
                                                  integrate_gradients,
                                                  this->subface_index);

  internal::FEFaceNormalEvaluationImpl<dim,
                                       fe_degree,
                                       n_components,
                                       VectorizedArray<Number>>::
    template interpolate<false, false>(
      *this->data, temp1, values_array, integrate_gradients, this->face_no);
}



template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number>
template <typename VectorType>
inline void
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  gather_evaluate(const VectorType &input_vector,
                  const bool        evaluate_values,
                  const bool        evaluate_gradients)
{
  const unsigned int side = this->face_no % 2;

  constexpr unsigned int static_dofs_per_face =
    fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
                     numbers::invalid_unsigned_int;
  const unsigned int dofs_per_face =
    fe_degree > -1 ? static_dofs_per_face :
                     Utilities::pow(this->data->fe_degree + 1, dim - 1);

  constexpr unsigned int stack_array_size_threshold = 100;

  VectorizedArray<Number>
    temp_data[static_dofs_per_face < stack_array_size_threshold ?
                n_components_ * 2 * dofs_per_face :
                1];
  VectorizedArray<Number> *__restrict temp1;
  if (static_dofs_per_face < stack_array_size_threshold)
    temp1 = &temp_data[0];
  else
    temp1 = this->scratch_data;

  internal::VectorReader<Number> reader;

  if (this->dof_info
          ->index_storage_variants[this->dof_access_index][this->cell] ==
        internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
          contiguous &&
      this->dof_info
          ->n_vectorization_lanes_filled[this->dof_access_index][this->cell] ==
        VectorizedArray<Number>::n_array_elements &&
      ((evaluate_gradients == false &&
        this->data->nodal_at_cell_boundaries == true) ||
       (this->data->element_type ==
          internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
        fe_degree > 1)))
    {
      const unsigned int *indices =
        &this->dof_info
           ->dof_indices_contiguous[this->dof_access_index]
                                   [this->cell *
                                    VectorizedArray<Number>::n_array_elements];
      if (evaluate_gradients == true &&
          this->data->element_type ==
            internal::MatrixFreeFunctions::tensor_symmetric_hermite)
        {
          // we know that the gradient weights for the Hermite case on the
          // right (side==1) are the negative from the value at the left
          // (side==0), so we only read out one of them.
          const VectorizedArray<Number> grad_weight0 =
            (side ? -1. : 1.) *
            this->data->shape_data_on_face[0][fe_degree + 1];
          const VectorizedArray<Number> grad_weight1 =
            (side ? -1. : 1.) *
            this->data->shape_data_on_face[0][fe_degree + 2];
          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
                          2 * dofs_per_face);

          const unsigned int *index_array =
            &this->data->face_to_cell_index_hermite(this->face_no, 0);
          for (unsigned int i = 0; i < dofs_per_face; ++i)
            {
              const unsigned int ind1 = index_array[2 * i];
              const unsigned int ind2 = index_array[2 * i + 1];
              for (unsigned int comp = 0; comp < n_components_; ++comp)
                {
                  reader.process_dof_gather(
                    indices,
                    input_vector,
                    ind1 + comp * static_dofs_per_component +
                      this->dof_info->component_dof_indices_offset
                        [this->active_fe_index][this->first_selected_component],
                    temp1[i + 2 * comp * dofs_per_face],
                    std::integral_constant<
                      bool,
                      std::is_same<typename VectorType::value_type,
                                   Number>::value>());
                  VectorizedArray<Number> grad;
                  reader.process_dof_gather(
                    indices,
                    input_vector,
                    ind2 + comp * static_dofs_per_component +
                      this->dof_info->component_dof_indices_offset
                        [this->active_fe_index][this->first_selected_component],
                    grad,
                    std::integral_constant<
                      bool,
                      std::is_same<typename VectorType::value_type,
                                   Number>::value>());
                  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
                    grad_weight0 * temp1[i + 2 * comp * dofs_per_face] +
                    grad_weight1 * grad;
                }
            }
        }
      else
        {
          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
                          dofs_per_face);
          const unsigned int *index_array =
            &this->data->face_to_cell_index_nodal(this->face_no, 0);
          for (unsigned int i = 0; i < dofs_per_face; ++i)
            for (unsigned int comp = 0; comp < n_components_; ++comp)
              {
                const unsigned int ind = index_array[i];
                reader.process_dof_gather(
                  indices,
                  input_vector,
                  ind + comp * static_dofs_per_component +
                    this->dof_info->component_dof_indices_offset
                      [this->active_fe_index][this->first_selected_component],
                  temp1[i + comp * 2 * dofs_per_face],
                  std::integral_constant<
                    bool,
                    std::is_same<typename VectorType::value_type,
                                 Number>::value>());
              }
        }
    }
  else
    {
      this->read_dof_values(input_vector);
      internal::FEFaceNormalEvaluationImpl<dim,
                                           fe_degree,
                                           n_components_,
                                           VectorizedArray<Number>>::
        template interpolate<true, false>(*this->data,
                                          this->values_dofs[0],
                                          temp1,
                                          evaluate_gradients,
                                          this->face_no);
    }

  if (fe_degree > -1 &&
      this->subface_index >= GeometryInfo<dim>::max_children_per_cell &&
      this->data->element_type <=
        internal::MatrixFreeFunctions::tensor_symmetric)
    internal::FEFaceEvaluationImpl<
      true,
      dim,
      fe_degree,
      n_q_points_1d,
      n_components_,
      VectorizedArray<Number>>::evaluate_in_face(*this->data,
                                                 temp1,
                                                 this->values_quad[0],
                                                 this->gradients_quad[0][0],
                                                 this->scratch_data +
                                                   2 * n_components_ *
                                                     dofs_per_face,
                                                 evaluate_values,
                                                 evaluate_gradients,
                                                 this->subface_index);
  else
    internal::FEFaceEvaluationImpl<
      false,
      dim,
      fe_degree,
      n_q_points_1d,
      n_components_,
      VectorizedArray<Number>>::evaluate_in_face(*this->data,
                                                 temp1,
                                                 this->values_quad[0],
                                                 this->gradients_quad[0][0],
                                                 this->scratch_data +
                                                   2 * n_components_ *
                                                     dofs_per_face,
                                                 evaluate_values,
                                                 evaluate_gradients,
                                                 this->subface_index);

  if (this->face_orientation)
    adjust_for_face_orientation(false, 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>
template <typename VectorType>
inline void
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  integrate_scatter(const bool  integrate_values,
                    const bool  integrate_gradients,
                    VectorType &destination)
{
  const unsigned int side = this->face_no % 2;
  const unsigned int dofs_per_face =
    fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
                     Utilities::pow(this->data->fe_degree + 1, dim - 1);

  constexpr unsigned int stack_array_size_threshold = 100;

  VectorizedArray<Number> temp_data[dofs_per_face < stack_array_size_threshold ?
                                      n_components_ * 2 * dofs_per_face :
                                      1];
  VectorizedArray<Number> *__restrict temp1;
  if (dofs_per_face < stack_array_size_threshold)
    temp1 = &temp_data[0];
  else
    temp1 = this->scratch_data;

  if (this->face_orientation)
    adjust_for_face_orientation(true, integrate_values, integrate_gradients);
  if (fe_degree > -1 &&
      this->subface_index >= GeometryInfo<dim>::max_children_per_cell &&
      this->data->element_type <=
        internal::MatrixFreeFunctions::tensor_symmetric)
    internal::FEFaceEvaluationImpl<
      true,
      dim,
      fe_degree,
      n_q_points_1d,
      n_components_,
      VectorizedArray<Number>>::integrate_in_face(*this->data,
                                                  temp1,
                                                  this->values_quad[0],
                                                  this->gradients_quad[0][0],
                                                  this->scratch_data +
                                                    2 * n_components_ *
                                                      dofs_per_face,
                                                  integrate_values,
                                                  integrate_gradients,
                                                  this->subface_index);
  else
    internal::FEFaceEvaluationImpl<
      false,
      dim,
      fe_degree,
      n_q_points_1d,
      n_components_,
      VectorizedArray<Number>>::integrate_in_face(*this->data,
                                                  temp1,
                                                  this->values_quad[0],
                                                  this->gradients_quad[0][0],
                                                  this->scratch_data +
                                                    2 * n_components_ *
                                                      dofs_per_face,
                                                  integrate_values,
                                                  integrate_gradients,
                                                  this->subface_index);

#  ifdef DEBUG
  this->dof_values_initialized = true;
#  endif

  internal::VectorDistributorLocalToGlobal<Number> writer;

  if (this->dof_info
          ->index_storage_variants[this->dof_access_index][this->cell] ==
        internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
          contiguous &&
      this->dof_info
          ->n_vectorization_lanes_filled[this->dof_access_index][this->cell] ==
        VectorizedArray<Number>::n_array_elements &&
      ((integrate_gradients == false &&
        this->data->nodal_at_cell_boundaries == true) ||
       (this->data->element_type ==
          internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
        fe_degree > 1)))
    {
      const unsigned int *indices =
        &this->dof_info
           ->dof_indices_contiguous[this->dof_access_index]
                                   [this->cell *
                                    VectorizedArray<Number>::n_array_elements];

      if (integrate_gradients == true &&
          this->data->element_type ==
            internal::MatrixFreeFunctions::tensor_symmetric_hermite)
        {
          // we know that the gradient weights for the Hermite case on the
          // right (side==1) are the negative from the value at the left
          // (side==0), so we only read out one of them.
          const VectorizedArray<Number> grad_weight0 =
            (side ? -1. : 1.) *
            this->data->shape_data_on_face[0][fe_degree + 1];
          const VectorizedArray<Number> grad_weight1 =
            (side ? -1. : 1.) *
            this->data->shape_data_on_face[0][fe_degree + 2];
          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
                          2 * dofs_per_face);
          const unsigned int *index_array =
            &this->data->face_to_cell_index_hermite(this->face_no, 0);
          for (unsigned int i = 0; i < dofs_per_face; ++i)
            {
              const unsigned int ind1 = index_array[2 * i];
              const unsigned int ind2 = index_array[2 * i + 1];
              for (unsigned int comp = 0; comp < n_components_; ++comp)
                {
                  VectorizedArray<Number> val =
                    temp1[i + 2 * comp * dofs_per_face] +
                    grad_weight0 *
                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
                  VectorizedArray<Number> grad =
                    grad_weight1 *
                    temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
                  writer.process_dof_gather(
                    indices,
                    destination,
                    comp * static_dofs_per_component + ind1 +
                      this->dof_info->component_dof_indices_offset
                        [this->active_fe_index][this->first_selected_component],
                    val,
                    std::integral_constant<
                      bool,
                      std::is_same<typename VectorType::value_type,
                                   Number>::value>());
                  writer.process_dof_gather(
                    indices,
                    destination,
                    comp * static_dofs_per_component + ind2 +
                      this->dof_info->component_dof_indices_offset
                        [this->active_fe_index][this->first_selected_component],
                    grad,
                    std::integral_constant<
                      bool,
                      std::is_same<typename VectorType::value_type,
                                   Number>::value>());
                }
            }
        }
      else
        {
          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
                          dofs_per_face);
          const unsigned int *index_array =
            &this->data->face_to_cell_index_nodal(this->face_no, 0);
          for (unsigned int i = 0; i < dofs_per_face; ++i)
            {
              const unsigned int ind = index_array[i];
              for (unsigned int comp = 0; comp < n_components_; ++comp)
                writer.process_dof_gather(
                  indices,
                  destination,
                  comp * static_dofs_per_component + ind +
                    this->dof_info->component_dof_indices_offset
                      [this->active_fe_index][this->first_selected_component],
                  temp1[i + 2 * comp * dofs_per_face],
                  std::integral_constant<
                    bool,
                    std::is_same<typename VectorType::value_type,
                                 Number>::value>());
            }
        }
    }
  else
    {
      internal::FEFaceNormalEvaluationImpl<dim,
                                           fe_degree,
                                           n_components_,
                                           VectorizedArray<Number>>::
        template interpolate<false, false>(*this->data,
                                           temp1,
                                           this->values_dofs[0],
                                           integrate_gradients,
                                           this->face_no);
      this->distribute_local_to_global(destination);
    }
}



template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components,
          typename Number>
inline void
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>::
  adjust_for_face_orientation(const bool integrate,
                              const bool values,
                              const bool gradients)
{
  VectorizedArray<Number> *tmp_values = this->scratch_data;
  const unsigned int *     orientations =
    &this->mapping_data->descriptor[this->active_fe_index]
       .face_orientations[this->face_orientation][0];
  for (unsigned int c = 0; c < n_components; ++c)
    {
      if (values == true)
        {
          if (integrate)
            for (unsigned int q = 0; q < n_q_points; ++q)
              tmp_values[orientations[q]] = this->values_quad[c][q];
          else
            for (unsigned int q = 0; q < n_q_points; ++q)
              tmp_values[q] = this->values_quad[c][orientations[q]];
          for (unsigned int q = 0; q < n_q_points; ++q)
            this->values_quad[c][q] = tmp_values[q];
        }
      if (gradients == true)
        for (unsigned int d = 0; d < dim; ++d)
          {
            if (integrate)
              for (unsigned int q = 0; q < n_q_points; ++q)
                tmp_values[orientations[q]] = this->gradients_quad[c][d][q];
            else
              for (unsigned int q = 0; q < n_q_points; ++q)
                tmp_values[q] = this->gradients_quad[c][d][orientations[q]];
            for (unsigned int q = 0; q < n_q_points; ++q)
              this->gradients_quad[c][d][q] = tmp_values[q];
          }
    }
}



template <int dim,
          int fe_degree,
          int n_q_points_1d,
          int n_components_,
          typename Number>
inline Point<dim, VectorizedArray<Number>>
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
  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