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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#ifndef dealii_matrix_free_fe_evaluation_h
#define dealii_matrix_free_fe_evaluation_h


#include <deal.II/base/array_view.h>
#include <deal.II/base/config.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/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/evaluation_kernels.h>
#include <deal.II/matrix_free/tensor_product_kernels.h>
#include <deal.II/matrix_free/evaluation_selector.h>

#include <deal.II/lac/vector_operation.h>

DEAL_II_NAMESPACE_OPEN



// forward declarations
namespace LinearAlgebra
{
  namespace distributed
  {
    template <typename> class Vector;
  }
}
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:
  typedef Number                            number_type;
  typedef Tensor<1,n_components_,VectorizedArray<Number> > value_type;
  typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
  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:
  typedef Number                            number_type;
  typedef Tensor<1,n_components_,VectorizedArray<Number> > value_type;
  typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
  static constexpr unsigned int dimension     = dim;
  static constexpr unsigned int n_components  = n_components_;
  typedef FEEvaluationBase<dim,n_components_, Number, is_face> BaseClass;

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:
  typedef Number                                 number_type;
  typedef VectorizedArray<Number>                value_type;
  typedef Tensor<1,dim,VectorizedArray<Number> > gradient_type;
  static constexpr unsigned int dimension          = dim;
  typedef FEEvaluationBase<dim,1,Number,is_face> BaseClass;

  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:
  typedef Number                            number_type;
  typedef Tensor<1,dim,VectorizedArray<Number> > value_type;
  typedef Tensor<2,dim,VectorizedArray<Number> > gradient_type;
  static constexpr unsigned int dimension     = dim;
  static constexpr unsigned int n_components  = dim;
  typedef FEEvaluationBase<dim,dim,Number,is_face> BaseClass;

  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:
  typedef Number                                 number_type;
  typedef VectorizedArray<Number>                value_type;
  typedef Tensor<1,1,VectorizedArray<Number> >   gradient_type;
  static constexpr unsigned int dimension      = 1;
  typedef FEEvaluationBase<1,1,Number,is_face>   BaseClass;

  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:
  typedef FEEvaluationAccess<dim,n_components_,Number,false> BaseClass;

  typedef Number                            number_type;

  typedef typename BaseClass::value_type    value_type;

  typedef typename BaseClass::gradient_type 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:
  typedef FEEvaluationAccess<dim,n_components_,Number,true> BaseClass;

  typedef Number                            number_type;

  typedef typename BaseClass::value_type    value_type;

  typedef typename BaseClass::gradient_type 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;
    };
  }
}


/*----------------------- 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>
  {
    typedef typename VectorType::BlockType BaseVectorType;

    static BaseVectorType *get_vector_component (VectorType &vec,
                                                 const unsigned int component)
    {
      AssertIndexRange (component, vec.n_blocks());
      return &vec.block(component);
    }
  };

  template <typename VectorType>
  struct BlockVectorSelector<VectorType,false>
  {
    typedef VectorType BaseVectorType;

    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>
  {
    typedef VectorType BaseVectorType;

    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>
  {
    typedef VectorType BaseVectorType;

    static BaseVectorType *get_vector_component (std::vector<VectorType *> &vec,
                                                 const unsigned int component)
    {
      AssertIndexRange (component, vec.size());
      return vec[component];
    }
  };
}



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)
                {
                  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]);
      }
  }
}



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;

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

  const 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;
  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;

  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