]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Simplify calls to FEEvaluation kernels
authorMartin Kronbichler <martin.kronbichler@it.uu.se>
Thu, 9 Dec 2021 17:13:08 +0000 (18:13 +0100)
committerMartin Kronbichler <martin.kronbichler@it.uu.se>
Sat, 11 Dec 2021 11:13:28 +0000 (12:13 +0100)
include/deal.II/fe/mapping_q.h
include/deal.II/fe/mapping_q_internal.h
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/evaluation_selector.h
include/deal.II/matrix_free/evaluation_template_factory.h
include/deal.II/matrix_free/evaluation_template_factory.templates.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/fe_evaluation_base_data.h
include/deal.II/matrix_free/mapping_info.templates.h
include/deal.II/matrix_free/mapping_info_storage.h
include/deal.II/matrix_free/mapping_info_storage.templates.h

index 306b257d5c64af50afa9ddbe56edfb70de791597..ddb9a633529ba004045484adeeae314e3157d723 100644 (file)
@@ -506,30 +506,6 @@ public:
      */
     mutable AlignedVector<VectorizedArrayType> scratch;
 
-    /**
-     * In case the quadrature rule given represents a tensor product
-     * the values at the mapped support points are stored in this object.
-     */
-    mutable AlignedVector<VectorizedArrayType> values_dofs;
-
-    /**
-     * In case the quadrature rule given represents a tensor product
-     * the values at the quadrature points are stored in this object.
-     */
-    mutable AlignedVector<VectorizedArrayType> values_quad;
-
-    /**
-     * In case the quadrature rule given represents a tensor product
-     * the gradients at the quadrature points are stored in this object.
-     */
-    mutable AlignedVector<VectorizedArrayType> gradients_quad;
-
-    /**
-     * In case the quadrature rule given represents a tensor product
-     * the hessians at the quadrature points are stored in this object.
-     */
-    mutable AlignedVector<VectorizedArrayType> hessians_quad;
-
     /**
      * Indicates whether the given Quadrature object is a tensor product.
      */
index 9df99a0404480d82cdaedaab865a84d599d6fd5e..3d40e5154d7780e3c96dbe6088df4689a85cf3b6 100644 (file)
@@ -35,6 +35,8 @@
 
 #include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/fe_evaluation_base_data.h>
+#include <deal.II/matrix_free/mapping_info_storage.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
 
@@ -1128,16 +1130,18 @@ namespace internal
           return;
         }
 
+      internal::MatrixFreeFunctions::
+        MappingInfoStorage<dim, dim, double, VectorizedArrayType>
+          temp_data;
+      temp_data.descriptor.resize(1);
+      temp_data.descriptor[0].n_q_points = n_q_points;
+      FEEvaluationBaseData<dim, double, false, VectorizedArrayType> eval(
+        {&data.shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0);
+
       // prepare arrays
       if (evaluation_flag != EvaluationFlags::nothing)
         {
-          data.values_dofs.resize(n_comp * n_shape_values);
-          data.values_quad.resize(n_comp * n_q_points);
-          data.gradients_quad.resize(n_comp * n_q_points * dim);
-          data.scratch.resize(2 * std::max(n_q_points, n_shape_values));
-
-          if (evaluation_flag & EvaluationFlags::hessians)
-            data.hessians_quad.resize(n_comp * n_q_points * n_hessians);
+          eval.set_data_pointers(&data.scratch, n_comp);
 
           const std::vector<unsigned int> &renumber_to_lexicographic =
             data.shape_info.lexicographic_numbering;
@@ -1146,20 +1150,14 @@ namespace internal
               {
                 const unsigned int in_comp  = d % n_lanes;
                 const unsigned int out_comp = d / n_lanes;
-                data.values_dofs[out_comp * n_shape_values + i][in_comp] =
+                eval
+                  .begin_dof_values()[out_comp * n_shape_values + i][in_comp] =
                   data.mapping_support_points[renumber_to_lexicographic[i]][d];
               }
 
           // do the actual tensorized evaluation
           internal::FEEvaluationFactory<dim, double, VectorizedArrayType>::
-            evaluate(n_comp,
-                     evaluation_flag,
-                     data.shape_info,
-                     data.values_dofs.begin(),
-                     data.values_quad.begin(),
-                     data.gradients_quad.begin(),
-                     data.hessians_quad.begin(),
-                     data.scratch.begin());
+            evaluate(n_comp, evaluation_flag, eval.begin_dof_values(), eval);
         }
 
       // do the postprocessing
@@ -1171,7 +1169,7 @@ namespace internal
                    in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
                    ++in_comp)
                 quadrature_points[i][out_comp * n_lanes + in_comp] =
-                  data.values_quad[out_comp * n_q_points + i][in_comp];
+                  eval.begin_values()[out_comp * n_q_points + i][in_comp];
         }
 
       if (evaluation_flag & EvaluationFlags::gradients)
@@ -1193,9 +1191,9 @@ namespace internal
                     const unsigned int new_point    = total_number % n_q_points;
                     data.contravariant[new_point][out_comp * n_lanes + in_comp]
                                       [new_comp] =
-                      data
-                        .gradients_quad[(out_comp * n_q_points + point) * dim +
-                                        j][in_comp];
+                      eval.begin_gradients()[(out_comp * n_q_points + point) *
+                                               dim +
+                                             j][in_comp];
                   }
         }
       if (update_flags & update_covariant_transformation)
@@ -1236,9 +1234,9 @@ namespace internal
                       dim == 2 ? desymmetrize_2d[new_hessian_comp][1] :
                                  desymmetrize_3d[new_hessian_comp][1];
                     const double value =
-                      data.hessians_quad[(out_comp * n_q_points + point) *
-                                           n_hessians +
-                                         j][in_comp];
+                      eval.begin_hessians()[(out_comp * n_q_points + point) *
+                                              n_hessians +
+                                            j][in_comp];
                     jacobian_grads[new_point][out_comp * n_lanes + in_comp]
                                   [new_hessian_comp_i][new_hessian_comp_j] =
                                     value;
index 81ff490f76840f0c9b2128870c1eba25fab69121..92001b04db438da77233dbb5f57cc4eaeb136bc3 100644 (file)
@@ -121,26 +121,21 @@ namespace internal
                                         n_q_points_1d,
                                         Number>;
 
+    using EvaluationData =
+      FEEvaluationBaseData<dim, typename Number::value_type, false, Number>;
+
     static void
-    evaluate(const unsigned int                            n_components,
-             const EvaluationFlags::EvaluationFlags        evaluation_flag,
-             const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-             const Number *                                values_dofs_actual,
-             Number *                                      values_quad,
-             Number *                                      gradients_quad,
-             Number *                                      hessians_quad,
-             Number *                                      scratch_data);
+    evaluate(const unsigned int                     n_components,
+             const EvaluationFlags::EvaluationFlags evaluation_flag,
+             const Number *                         values_dofs_actual,
+             EvaluationData &                       eval);
 
     static void
-    integrate(const unsigned int                            n_components,
-              const EvaluationFlags::EvaluationFlags        integration_flag,
-              const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-              Number *                                      values_dofs_actual,
-              Number *                                      values_quad,
-              Number *                                      gradients_quad,
-              Number *                                      hessians_quad,
-              Number *                                      scratch_data,
-              const bool add_into_values_array);
+    integrate(const unsigned int                     n_components,
+              const EvaluationFlags::EvaluationFlags integration_flag,
+              Number *                               values_dofs_actual,
+              EvaluationData &                       eval,
+              const bool                             add_into_values_array);
 
     static Eval
     create_evaluator_tensor_product(
@@ -175,26 +170,21 @@ namespace internal
                           n_q_points_1d,
                           Number>
   {
+    using EvaluationData =
+      FEEvaluationBaseData<dim, typename Number::value_type, false, Number>;
+
     static void
-    evaluate(const unsigned int                            n_components,
-             const EvaluationFlags::EvaluationFlags        evaluation_flag,
-             const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-             const Number *                                values_dofs_actual,
-             Number *                                      values_quad,
-             Number *                                      gradients_quad,
-             Number *                                      hessians_quad,
-             Number *                                      scratch_data);
+    evaluate(const unsigned int                     n_components,
+             const EvaluationFlags::EvaluationFlags evaluation_flag,
+             const Number *                         values_dofs_actual,
+             EvaluationData &                       eval);
 
     static void
-    integrate(const unsigned int                            n_components,
-              const EvaluationFlags::EvaluationFlags        integration_flag,
-              const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-              Number *                                      values_dofs_actual,
-              Number *                                      values_quad,
-              Number *                                      gradients_quad,
-              Number *                                      hessians_quad,
-              Number *                                      scratch_data,
-              const bool add_into_values_array);
+    integrate(const unsigned int                     n_components,
+              const EvaluationFlags::EvaluationFlags integration_flag,
+              Number *                               values_dofs_actual,
+              EvaluationData &                       eval,
+              const bool                             add_into_values_array);
   };
 
 
@@ -206,14 +196,10 @@ namespace internal
             typename Number>
   inline void
   FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, Number>::evaluate(
-    const unsigned int                            n_components,
-    const EvaluationFlags::EvaluationFlags        evaluation_flag,
-    const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-    const Number *                                values_dofs_actual,
-    Number *                                      values_quad,
-    Number *                                      gradients_quad,
-    Number *                                      hessians_quad,
-    Number *                                      scratch_data)
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags evaluation_flag,
+    const Number *                         values_dofs_actual,
+    EvaluationData &                       eval)
   {
     if (evaluation_flag == EvaluationFlags::nothing)
       return;
@@ -221,11 +207,13 @@ namespace internal
     std::array<const MatrixFreeFunctions::UnivariateShapeData<Number> *, 3>
       univariate_shape_data;
 
-    univariate_shape_data.fill(&shape_info.data[0]);
+    const auto &shape_data = eval.get_shape_info().data;
 
-    if (shape_info.data.size() == dim)
+    univariate_shape_data.fill(&shape_data[0]);
+
+    if (shape_data.size() == dim)
       for (int i = 1; i < dim; ++i)
-        univariate_shape_data[i] = &shape_info.data[i];
+        univariate_shape_data[i] = &shape_data[i];
 
     Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
     Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
@@ -237,34 +225,36 @@ namespace internal
         (Eval::n_rows_of_product > Eval::n_columns_of_product ?
            Eval::n_rows_of_product :
            Eval::n_columns_of_product);
-    Number *temp1 = scratch_data;
+    Number *temp1 = eval.get_scratch_data().begin();
     Number *temp2;
     if (temp_size == 0)
       {
-        temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
-                                   shape_info.data.front().fe_degree + 1),
-                                 Utilities::fixed_power<dim>(
-                                   shape_info.data.front().n_q_points_1d));
+        temp2 =
+          temp1 +
+          std::max(Utilities::fixed_power<dim>(shape_data[0].fe_degree + 1),
+                   Utilities::fixed_power<dim>(shape_data[0].n_q_points_1d));
       }
     else
       {
         temp2 = temp1 + temp_size;
       }
 
-    const unsigned int n_q_points =
-      temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
-    const unsigned int dofs_per_comp =
+    const std::size_t n_q_points = temp_size == 0 ?
+                                     eval.get_shape_info().n_q_points :
+                                     Eval::n_columns_of_product;
+    const std::size_t dofs_per_comp =
       (type == MatrixFreeFunctions::truncated_tensor) ?
-        Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
-        shape_info.dofs_per_component_on_cell;
+        Utilities::pow(shape_data[0].fe_degree + 1, dim) :
+        eval.get_shape_info().dofs_per_component_on_cell;
     const Number *values_dofs = values_dofs_actual;
     if (type == MatrixFreeFunctions::truncated_tensor)
       {
+        const std::size_t stored_dofs_per_comp =
+          eval.get_shape_info().dofs_per_component_on_cell;
         Number *values_dofs_tmp =
-          scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
-                                       shape_info.n_q_points));
+          temp1 + 2 * (std::max(stored_dofs_per_comp, n_q_points));
         const int degree =
-          fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
+          fe_degree != -1 ? fe_degree : shape_data[0].fe_degree;
         for (unsigned int c = 0; c < n_components; ++c)
           for (int i = 0, count_p = 0, count_q = 0;
                i < (dim > 2 ? degree + 1 : 1);
@@ -275,8 +265,7 @@ namespace internal
                   for (int k = 0; k < degree + 1 - j - i;
                        ++k, ++count_p, ++count_q)
                     values_dofs_tmp[c * dofs_per_comp + count_q] =
-                      values_dofs_actual
-                        [c * shape_info.dofs_per_component_on_cell + count_p];
+                      values_dofs_actual[c * stored_dofs_per_comp + count_p];
                   for (int k = degree + 1 - j - i; k < degree + 1;
                        ++k, ++count_q)
                     values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
@@ -288,6 +277,10 @@ namespace internal
         values_dofs = values_dofs_tmp;
       }
 
+    Number *values_quad    = eval.begin_values();
+    Number *gradients_quad = eval.begin_gradients();
+    Number *hessians_quad  = eval.begin_hessians();
+
     switch (dim)
       {
         case 1:
@@ -460,10 +453,10 @@ namespace internal
       {
         values_quad -= n_components * n_q_points;
         values_dofs -= n_components * dofs_per_comp;
-        for (unsigned int c = 0; c < n_components; ++c)
-          for (unsigned int q = 0; q < shape_info.n_q_points; ++q)
-            values_quad[c * shape_info.n_q_points + q] +=
-              values_dofs[(c + 1) * shape_info.dofs_per_component_on_cell - 1];
+        for (std::size_t c = 0; c < n_components; ++c)
+          for (std::size_t q = 0; q < n_q_points; ++q)
+            values_quad[c * n_q_points + q] +=
+              values_dofs[(c + 1) * dofs_per_comp - 1];
       }
   }
 
@@ -476,24 +469,21 @@ namespace internal
             typename Number>
   inline void
   FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, Number>::integrate(
-    const unsigned int                            n_components,
-    const EvaluationFlags::EvaluationFlags        integration_flag,
-    const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-    Number *                                      values_dofs_actual,
-    Number *                                      values_quad,
-    Number *                                      gradients_quad,
-    Number *                                      hessians_quad,
-    Number *                                      scratch_data,
-    const bool                                    add_into_values_array)
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags integration_flag,
+    Number *                               values_dofs_actual,
+    EvaluationData &                       eval,
+    const bool                             add_into_values_array)
   {
     std::array<const MatrixFreeFunctions::UnivariateShapeData<Number> *, 3>
       univariate_shape_data;
 
-    univariate_shape_data.fill(&shape_info.data[0]);
+    const auto &shape_data = eval.get_shape_info().data;
+    univariate_shape_data.fill(&shape_data[0]);
 
-    if (shape_info.data.size() == dim)
+    if (shape_data.size() == dim)
       for (int i = 1; i < dim; ++i)
-        univariate_shape_data[i] = &shape_info.data[i];
+        univariate_shape_data[i] = &shape_data[i];
 
     Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
     Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
@@ -505,33 +495,38 @@ namespace internal
         (Eval::n_rows_of_product > Eval::n_columns_of_product ?
            Eval::n_rows_of_product :
            Eval::n_columns_of_product);
-    Number *temp1 = scratch_data;
+    Number *temp1 = eval.get_scratch_data().begin();
     Number *temp2;
     if (temp_size == 0)
       {
-        temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
-                                   shape_info.data.front().fe_degree + 1),
-                                 Utilities::fixed_power<dim>(
-                                   shape_info.data.front().n_q_points_1d));
+        temp2 =
+          temp1 +
+          std::max(Utilities::fixed_power<dim>(shape_data[0].fe_degree + 1),
+                   Utilities::fixed_power<dim>(shape_data[0].n_q_points_1d));
       }
     else
       {
         temp2 = temp1 + temp_size;
       }
 
-    const unsigned int n_q_points =
-      temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
+    const unsigned int n_q_points = temp_size == 0 ?
+                                      eval.get_shape_info().n_q_points :
+                                      Eval::n_columns_of_product;
     const unsigned int dofs_per_comp =
       (type == MatrixFreeFunctions::truncated_tensor) ?
-        Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
-        shape_info.dofs_per_component_on_cell;
+        Utilities::fixed_power<dim>(shape_data[0].fe_degree + 1) :
+        eval.get_shape_info().dofs_per_component_on_cell;
     // expand dof_values to tensor product for truncated tensor products
     Number *values_dofs =
       (type == MatrixFreeFunctions::truncated_tensor) ?
-        scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
-                                     shape_info.n_q_points)) :
+        temp1 + 2 * (std::max(eval.get_shape_info().dofs_per_component_on_cell,
+                              n_q_points)) :
         values_dofs_actual;
 
+    Number *values_quad    = eval.begin_values();
+    Number *gradients_quad = eval.begin_gradients();
+    Number *hessians_quad  = eval.begin_hessians();
+
     switch (dim)
       {
         case 1:
@@ -734,14 +729,13 @@ namespace internal
     // case FE_Q_DG0: add values, gradients and second derivatives are zero
     if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
       {
-        values_dofs -= n_components * dofs_per_comp -
-                       shape_info.dofs_per_component_on_cell + 1;
+        values_dofs -= n_components * dofs_per_comp - dofs_per_comp + 1;
         values_quad -= n_components * n_q_points;
         if (integration_flag & EvaluationFlags::values)
           for (unsigned int c = 0; c < n_components; ++c)
             {
               values_dofs[0] = values_quad[0];
-              for (unsigned int q = 1; q < shape_info.n_q_points; ++q)
+              for (unsigned int q = 1; q < n_q_points; ++q)
                 values_dofs[0] += values_quad[q];
               values_dofs += dofs_per_comp;
               values_quad += n_q_points;
@@ -749,16 +743,18 @@ namespace internal
         else
           {
             for (unsigned int c = 0; c < n_components; ++c)
-              values_dofs[c * shape_info.dofs_per_component_on_cell] = Number();
-            values_dofs += n_components * shape_info.dofs_per_component_on_cell;
+              values_dofs[c * dofs_per_comp] = Number();
+            values_dofs += n_components * dofs_per_comp;
           }
       }
 
     if (type == MatrixFreeFunctions::truncated_tensor)
       {
+        const std::size_t stored_dofs_per_comp =
+          eval.get_shape_info().dofs_per_component_on_cell;
         values_dofs -= dofs_per_comp * n_components;
         const int degree =
-          fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
+          fe_degree != -1 ? fe_degree : shape_data[0].fe_degree;
         for (unsigned int c = 0; c < n_components; ++c)
           for (int i = 0, count_p = 0, count_q = 0;
                i < (dim > 2 ? degree + 1 : 1);
@@ -768,9 +764,7 @@ namespace internal
                 {
                   for (int k = 0; k < degree + 1 - j - i;
                        ++k, ++count_p, ++count_q)
-                    values_dofs_actual[c *
-                                         shape_info.dofs_per_component_on_cell +
-                                       count_p] =
+                    values_dofs_actual[c * stored_dofs_per_comp + count_p] =
                       values_dofs[c * dofs_per_comp + count_q];
                   count_q += j + i;
                 }
@@ -783,32 +777,29 @@ namespace internal
 
   template <int dim, int fe_degree, int n_q_points_1d, typename Number>
   inline void
-  FEEvaluationImpl<
-    MatrixFreeFunctions::tensor_none,
-    dim,
-    fe_degree,
-    n_q_points_1d,
-    Number>::evaluate(const unsigned int                     n_components,
-                      const EvaluationFlags::EvaluationFlags evaluation_flag,
-                      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-                      const Number *values_dofs_actual,
-                      Number *      values_quad,
-                      Number *      gradients_quad,
-                      Number *      hessians_quad,
-                      Number *      scratch_data)
+  FEEvaluationImpl<MatrixFreeFunctions::tensor_none,
+                   dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   Number>::evaluate(const unsigned int n_components,
+                                     const EvaluationFlags::EvaluationFlags
+                                                     evaluation_flag,
+                                     const Number *  values_dofs_actual,
+                                     EvaluationData &eval)
   {
-    (void)scratch_data;
+    const unsigned int n_dofs =
+      eval.get_shape_info().dofs_per_component_on_cell;
+    const unsigned int n_q_points = eval.get_shape_info().n_q_points;
 
-    const unsigned int n_dofs     = shape_info.dofs_per_component_on_cell;
-    const unsigned int n_q_points = shape_info.n_q_points;
+    const auto &shape_data = eval.get_shape_info().data;
 
     using Eval =
       EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
 
     if (evaluation_flag & EvaluationFlags::values)
       {
-        const auto shape_values = shape_info.data.front().shape_values.data();
-        auto       values_quad_ptr        = values_quad;
+        const auto shape_values           = shape_data[0].shape_values.data();
+        auto       values_quad_ptr        = eval.begin_values();
         auto       values_dofs_actual_ptr = values_dofs_actual;
 
         Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
@@ -824,10 +815,9 @@ namespace internal
 
     if (evaluation_flag & EvaluationFlags::gradients)
       {
-        const auto shape_gradients =
-          shape_info.data.front().shape_gradients.data();
-        auto gradients_quad_ptr     = gradients_quad;
-        auto values_dofs_actual_ptr = values_dofs_actual;
+        const auto shape_gradients    = shape_data[0].shape_gradients.data();
+        auto       gradients_quad_ptr = eval.begin_gradients();
+        auto       values_dofs_actual_ptr = values_dofs_actual;
 
         for (unsigned int c = 0; c < n_components; ++c)
           {
@@ -849,47 +839,41 @@ namespace internal
       }
 
     if (evaluation_flag & EvaluationFlags::hessians)
-      {
-        Assert(false, ExcNotImplemented());
-        (void)hessians_quad;
-      }
+      Assert(false, ExcNotImplemented());
   }
 
 
 
   template <int dim, int fe_degree, int n_q_points_1d, typename Number>
   inline void
-  FEEvaluationImpl<
-    MatrixFreeFunctions::tensor_none,
-    dim,
-    fe_degree,
-    n_q_points_1d,
-    Number>::integrate(const unsigned int                     n_components,
-                       const EvaluationFlags::EvaluationFlags integration_flag,
-                       const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-                       Number *   values_dofs_actual,
-                       Number *   values_quad,
-                       Number *   gradients_quad,
-                       Number *   hessians_quad,
-                       Number *   scratch_data,
-                       const bool add_into_values_array)
+  FEEvaluationImpl<MatrixFreeFunctions::tensor_none,
+                   dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   Number>::integrate(const unsigned int n_components,
+                                      const EvaluationFlags::EvaluationFlags
+                                                      integration_flag,
+                                      Number *        values_dofs_actual,
+                                      EvaluationData &eval,
+                                      const bool      add_into_values_array)
   {
     // TODO: implement hessians
-    (void)hessians_quad;
     AssertThrow(!(integration_flag & EvaluationFlags::hessians),
                 ExcNotImplemented());
-    (void)scratch_data;
 
-    const unsigned int n_dofs     = shape_info.dofs_per_component_on_cell;
-    const unsigned int n_q_points = shape_info.n_q_points;
+    const unsigned int n_dofs =
+      eval.get_shape_info().dofs_per_component_on_cell;
+    const unsigned int n_q_points = eval.get_shape_info().n_q_points;
+
+    const auto &shape_data = eval.get_shape_info().data;
 
     using Eval =
       EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
 
     if (integration_flag & EvaluationFlags::values)
       {
-        const auto shape_values = shape_info.data.front().shape_values.data();
-        auto       values_quad_ptr        = values_quad;
+        const auto shape_values           = shape_data[0].shape_values.data();
+        auto       values_quad_ptr        = eval.begin_values();
         auto       values_dofs_actual_ptr = values_dofs_actual;
 
         Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
@@ -909,10 +893,9 @@ namespace internal
 
     if (integration_flag & EvaluationFlags::gradients)
       {
-        const auto shape_gradients =
-          shape_info.data.front().shape_gradients.data();
-        auto gradients_quad_ptr     = gradients_quad;
-        auto values_dofs_actual_ptr = values_dofs_actual;
+        const auto shape_gradients    = shape_data[0].shape_gradients.data();
+        auto       gradients_quad_ptr = eval.begin_gradients();
+        auto       values_dofs_actual_ptr = values_dofs_actual;
 
         for (unsigned int c = 0; c < n_components; ++c)
           {
@@ -1357,26 +1340,36 @@ namespace internal
   template <int dim, int fe_degree, typename Number>
   struct FEEvaluationImplCollocation
   {
+    using EvaluationData =
+      FEEvaluationBaseData<dim, typename Number::value_type, false, Number>;
+
+    static void
+    evaluate(const unsigned int                     n_components,
+             const EvaluationFlags::EvaluationFlags evaluation_flag,
+             const Number *                         values_dofs,
+             EvaluationData &                       eval);
+
+    static void
+    do_evaluate(const MatrixFreeFunctions::UnivariateShapeData<Number> &shape,
+                const EvaluationFlags::EvaluationFlags evaluation_flag,
+                const Number *                         values_dofs,
+                Number *                               gradients_quad,
+                Number *                               hessians_quad);
+
     static void
-    evaluate(const unsigned int                            n_components,
-             const EvaluationFlags::EvaluationFlags        evaluation_flag,
-             const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-             const Number *                                values_dofs,
-             Number *                                      values_quad,
-             Number *                                      gradients_quad,
-             Number *                                      hessians_quad,
-             Number *                                      scratch_data);
+    integrate(const unsigned int                     n_components,
+              const EvaluationFlags::EvaluationFlags integration_flag,
+              Number *                               values_dofs,
+              EvaluationData &                       eval,
+              const bool                             add_into_values_array);
 
     static void
-    integrate(const unsigned int                            n_components,
-              const EvaluationFlags::EvaluationFlags        integration_flag,
-              const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-              Number *                                      values_dofs,
-              Number *                                      values_quad,
-              Number *                                      gradients_quad,
-              Number *                                      hessians_quad,
-              Number *                                      scratch_data,
-              const bool add_into_values_array);
+    do_integrate(const MatrixFreeFunctions::UnivariateShapeData<Number> &shape,
+                 const EvaluationFlags::EvaluationFlags integration_flag,
+                 Number *                               values_dofs,
+                 Number *                               gradients_quad,
+                 const Number *                         hessians_quad,
+                 const bool                             add_into_values_array);
   };
 
 
@@ -1384,76 +1377,86 @@ namespace internal
   template <int dim, int fe_degree, typename Number>
   inline void
   FEEvaluationImplCollocation<dim, fe_degree, Number>::evaluate(
-    const unsigned int                            n_components,
-    const EvaluationFlags::EvaluationFlags        evaluation_flag,
-    const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-    const Number *                                values_dofs,
-    Number *                                      values_quad,
-    Number *                                      gradients_quad,
-    Number *                                      hessians_quad,
-    Number *)
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags evaluation_flag,
+    const Number *                         values_dofs,
+    EvaluationData &                       eval)
+  {
+    constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
+
+    for (unsigned int c = 0; c < n_components; ++c)
+      {
+        if (evaluation_flag & EvaluationFlags::values)
+          for (unsigned int i = 0; i < n_points; ++i)
+            eval.begin_values()[n_points * c + i] =
+              values_dofs[n_points * c + i];
+
+        do_evaluate(eval.get_shape_info().data[0],
+                    evaluation_flag,
+                    values_dofs + c * n_points,
+                    eval.begin_gradients() + c * dim * n_points,
+                    eval.begin_hessians() + c * dim * (dim + 1) / 2 * n_points);
+      }
+  }
+
+
+
+  template <int dim, int fe_degree, typename Number>
+  inline void
+  FEEvaluationImplCollocation<dim, fe_degree, Number>::do_evaluate(
+    const MatrixFreeFunctions::UnivariateShapeData<Number> &shape,
+    const EvaluationFlags::EvaluationFlags                  evaluation_flag,
+    const Number *                                          values_dofs,
+    Number *                                                gradients_quad,
+    Number *                                                hessians_quad)
   {
-    AssertDimension(
-      shape_info.data.front().shape_gradients_collocation_eo.size(),
-      (fe_degree + 2) / 2 * (fe_degree + 1));
+    AssertDimension(shape.shape_gradients_collocation_eo.size(),
+                    (fe_degree + 2) / 2 * (fe_degree + 1));
+    constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
 
     EvaluatorTensorProduct<evaluate_evenodd,
                            dim,
                            fe_degree + 1,
                            fe_degree + 1,
                            Number>
-                           eval(AlignedVector<Number>(),
-           shape_info.data.front().shape_gradients_collocation_eo,
-           shape_info.data.front().shape_hessians_collocation_eo);
-    constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
-
-    for (unsigned int c = 0; c < n_components; ++c)
+      eval(AlignedVector<Number>(),
+           shape.shape_gradients_collocation_eo,
+           shape.shape_hessians_collocation_eo);
+    if (evaluation_flag &
+        (EvaluationFlags::gradients | EvaluationFlags::hessians))
       {
-        if (evaluation_flag & EvaluationFlags::values)
-          for (unsigned int i = 0; i < n_q_points; ++i)
-            values_quad[i] = values_dofs[i];
-        if (evaluation_flag &
-            (EvaluationFlags::gradients | EvaluationFlags::hessians))
+        eval.template gradients<0, true, false>(values_dofs, gradients_quad);
+        if (dim > 1)
+          eval.template gradients<1, true, false>(values_dofs,
+                                                  gradients_quad + n_points);
+        if (dim > 2)
+          eval.template gradients<2, true, false>(values_dofs,
+                                                  gradients_quad +
+                                                    2 * n_points);
+      }
+    if (evaluation_flag & EvaluationFlags::hessians)
+      {
+        eval.template hessians<0, true, false>(values_dofs, hessians_quad);
+        if (dim > 1)
           {
-            eval.template gradients<0, true, false>(values_dofs,
-                                                    gradients_quad);
-            if (dim > 1)
-              eval.template gradients<1, true, false>(values_dofs,
-                                                      gradients_quad +
-                                                        n_q_points);
-            if (dim > 2)
-              eval.template gradients<2, true, false>(values_dofs,
-                                                      gradients_quad +
-                                                        2 * n_q_points);
+            eval.template gradients<1, true, false>(gradients_quad,
+                                                    hessians_quad +
+                                                      dim * n_points);
+            eval.template hessians<1, true, false>(values_dofs,
+                                                   hessians_quad + n_points);
           }
-        if (evaluation_flag & EvaluationFlags::hessians)
+        if (dim > 2)
           {
-            eval.template hessians<0, true, false>(values_dofs, hessians_quad);
-            if (dim > 1)
-              {
-                eval.template gradients<1, true, false>(gradients_quad,
-                                                        hessians_quad +
-                                                          dim * n_q_points);
-                eval.template hessians<1, true, false>(values_dofs,
-                                                       hessians_quad +
-                                                         n_q_points);
-              }
-            if (dim > 2)
-              {
-                eval.template gradients<2, true, false>(gradients_quad,
-                                                        hessians_quad +
-                                                          4 * n_q_points);
-                eval.template gradients<2, true, false>(
-                  gradients_quad + n_q_points, hessians_quad + 5 * n_q_points);
-                eval.template hessians<2, true, false>(values_dofs,
-                                                       hessians_quad +
-                                                         2 * n_q_points);
-              }
-            hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
+            eval.template gradients<2, true, false>(gradients_quad,
+                                                    hessians_quad +
+                                                      4 * n_points);
+            eval.template gradients<2, true, false>(gradients_quad + n_points,
+                                                    hessians_quad +
+                                                      5 * n_points);
+            eval.template hessians<2, true, false>(values_dofs,
+                                                   hessians_quad +
+                                                     2 * n_points);
           }
-        gradients_quad += dim * n_q_points;
-        values_quad += n_q_points;
-        values_dofs += n_q_points;
       }
   }
 
@@ -1462,117 +1465,140 @@ namespace internal
   template <int dim, int fe_degree, typename Number>
   inline void
   FEEvaluationImplCollocation<dim, fe_degree, Number>::integrate(
-    const unsigned int                            n_components,
-    const EvaluationFlags::EvaluationFlags        integration_flag,
-    const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-    Number *                                      values_dofs,
-    Number *                                      values_quad,
-    Number *                                      gradients_quad,
-    Number *                                      hessians_quad,
-    Number *                                      scratch_data,
-    const bool                                    add_into_values_array)
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags integration_flag,
+    Number *                               values_dofs,
+    EvaluationData &                       eval,
+    const bool                             add_into_values_array)
+  {
+    constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
+
+    for (unsigned int c = 0; c < n_components; ++c)
+      {
+        if (integration_flag & EvaluationFlags::values)
+          {
+            if (add_into_values_array)
+              for (unsigned int i = 0; i < n_points; ++i)
+                values_dofs[n_points * c + i] +=
+                  eval.begin_values()[n_points * c + i];
+            else
+              for (unsigned int i = 0; i < n_points; ++i)
+                values_dofs[n_points * c + i] =
+                  eval.begin_values()[n_points * c + i];
+          }
+
+        do_integrate(eval.get_shape_info().data[0],
+                     integration_flag,
+                     values_dofs + c * n_points,
+                     eval.begin_gradients() + c * dim * n_points,
+                     eval.begin_hessians() + c * dim * (dim + 1) / 2 * n_points,
+                     add_into_values_array ||
+                       (integration_flag & EvaluationFlags::values));
+      }
+  }
+
+
+
+  template <int dim, int fe_degree, typename Number>
+  inline void
+  FEEvaluationImplCollocation<dim, fe_degree, Number>::do_integrate(
+    const MatrixFreeFunctions::UnivariateShapeData<Number> &shape,
+    const EvaluationFlags::EvaluationFlags                  integration_flag,
+    Number *                                                values_dofs,
+    Number *                                                gradients_quad,
+    const Number *                                          hessians_quad,
+    const bool add_into_values_array)
   {
-    AssertDimension(
-      shape_info.data.front().shape_gradients_collocation_eo.size(),
-      (fe_degree + 2) / 2 * (fe_degree + 1));
+    AssertDimension(shape.shape_gradients_collocation_eo.size(),
+                    (fe_degree + 2) / 2 * (fe_degree + 1));
 
     EvaluatorTensorProduct<evaluate_evenodd,
                            dim,
                            fe_degree + 1,
                            fe_degree + 1,
                            Number>
-                           eval(AlignedVector<Number>(),
-           shape_info.data.front().shape_gradients_collocation_eo,
-           shape_info.data.front().shape_hessians_collocation_eo);
-    constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
-    constexpr unsigned int hdim       = (dim * (dim + 1)) / 2;
+                          eval(AlignedVector<Number>(),
+           shape.shape_gradients_collocation_eo,
+           shape.shape_hessians_collocation_eo);
+    constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
 
-    for (unsigned int c = 0; c < n_components; ++c)
+    if (integration_flag & EvaluationFlags::hessians)
       {
-        if (integration_flag & EvaluationFlags::values)
+        // diagonal
+        // grad xx
+        if (add_into_values_array == true)
+          eval.template hessians<0, false, true>(hessians_quad, values_dofs);
+        else
+          eval.template hessians<0, false, false>(hessians_quad, values_dofs);
+        // grad yy
+        if (dim > 1)
+          eval.template hessians<1, false, true>(hessians_quad + n_points,
+                                                 values_dofs);
+        // grad zz
+        if (dim > 2)
+          eval.template hessians<2, false, true>(hessians_quad + 2 * n_points,
+                                                 values_dofs);
+        // off-diagonal
+        if (dim == 2)
           {
-            if (add_into_values_array == false)
-              for (unsigned int i = 0; i < n_q_points; ++i)
-                values_dofs[i] = values_quad[i];
+            // grad xy, queue into gradient
+            if (integration_flag & EvaluationFlags::gradients)
+              eval.template gradients<1, false, true>(hessians_quad +
+                                                        2 * n_points,
+                                                      gradients_quad);
             else
-              for (unsigned int i = 0; i < n_q_points; ++i)
-                values_dofs[i] += values_quad[i];
+              eval.template gradients<1, false, false>(hessians_quad +
+                                                         2 * n_points,
+                                                       gradients_quad);
           }
-        if (integration_flag & EvaluationFlags::gradients)
+        if (dim == 3)
           {
-            if (integration_flag & EvaluationFlags::values ||
-                add_into_values_array == true)
-              eval.template gradients<0, false, true>(gradients_quad,
-                                                      values_dofs);
+            // grad xy, queue into gradient
+            if (integration_flag & EvaluationFlags::gradients)
+              eval.template gradients<1, false, true>(hessians_quad +
+                                                        3 * n_points,
+                                                      gradients_quad);
             else
-              eval.template gradients<0, false, false>(gradients_quad,
-                                                       values_dofs);
-            if (dim > 1)
-              eval.template gradients<1, false, true>(gradients_quad +
-                                                        n_q_points,
-                                                      values_dofs);
-            if (dim > 2)
-              eval.template gradients<2, false, true>(gradients_quad +
-                                                        2 * n_q_points,
-                                                      values_dofs);
-          }
-        if (integration_flag & EvaluationFlags::hessians)
-          {
-            // diagonal
-            // grad xx
-            if (integration_flag & EvaluationFlags::values ||
-                integration_flag & EvaluationFlags::gradients ||
-                add_into_values_array == true)
-              eval.template hessians<0, false, true>(hessians_quad,
-                                                     values_dofs);
+              eval.template gradients<1, false, false>(hessians_quad +
+                                                         3 * n_points,
+                                                       gradients_quad);
+
+            // grad xz
+            eval.template gradients<2, false, true>(hessians_quad +
+                                                      4 * n_points,
+                                                    gradients_quad);
+
+            // grad yz
+            if (integration_flag & EvaluationFlags::gradients)
+              eval.template gradients<2, false, true>(
+                hessians_quad + 5 * n_points, gradients_quad + n_points);
             else
-              eval.template hessians<0, false, false>(hessians_quad,
-                                                      values_dofs);
-            // grad yy
-            if (dim > 1)
-              eval.template hessians<1, false, true>(hessians_quad + n_q_points,
-                                                     values_dofs);
-            // grad zz
-            if (dim > 2)
-              eval.template hessians<2, false, true>(hessians_quad +
-                                                       2 * n_q_points,
-                                                     values_dofs);
-            // off-diagonal
-            if (dim == 2)
-              {
-                // grad xy
-                eval.template gradients<0, false, false>(hessians_quad +
-                                                           2 * n_q_points,
-                                                         scratch_data);
-                eval.template gradients<1, false, true>(scratch_data,
-                                                        values_dofs);
-              }
-            if (dim == 3)
-              {
-                // grad xy
-                eval.template gradients<0, false, false>(hessians_quad +
-                                                           3 * n_q_points,
-                                                         scratch_data);
-                eval.template gradients<1, false, true>(scratch_data,
-                                                        values_dofs);
-                // grad xz
-                eval.template gradients<0, false, false>(hessians_quad +
-                                                           4 * n_q_points,
-                                                         scratch_data);
-                eval.template gradients<2, false, true>(scratch_data,
-                                                        values_dofs);
-                // grad yz
-                eval.template gradients<1, false, false>(hessians_quad +
-                                                           5 * n_q_points,
-                                                         scratch_data);
-                eval.template gradients<2, false, true>(scratch_data,
-                                                        values_dofs);
-              }
-            hessians_quad += hdim * n_q_points;
+              eval.template gradients<2, false, false>(
+                hessians_quad + 5 * n_points, gradients_quad + n_points);
           }
-        gradients_quad += dim * n_q_points;
-        values_quad += n_q_points;
-        values_dofs += n_q_points;
+
+        // if we did not integrate gradients, set the last slot to zero
+        // which was not touched before, in order to avoid the if
+        // statement in the gradients loop below
+        if (!(integration_flag & EvaluationFlags::gradients))
+          for (unsigned int q = 0; q < n_points; ++q)
+            gradients_quad[(dim - 1) * n_points + q] = Number();
+      }
+
+    if (integration_flag &
+        (EvaluationFlags::gradients | EvaluationFlags::hessians))
+      {
+        if (add_into_values_array ||
+            integration_flag & EvaluationFlags::hessians)
+          eval.template gradients<0, false, true>(gradients_quad, values_dofs);
+        else
+          eval.template gradients<0, false, false>(gradients_quad, values_dofs);
+        if (dim > 1)
+          eval.template gradients<1, false, true>(gradients_quad + n_points,
+                                                  values_dofs);
+        if (dim > 2)
+          eval.template gradients<2, false, true>(gradients_quad + 2 * n_points,
+                                                  values_dofs);
       }
   }
 
@@ -1591,26 +1617,21 @@ namespace internal
   template <int dim, int fe_degree, int n_q_points_1d, typename Number>
   struct FEEvaluationImplTransformToCollocation
   {
+    using EvaluationData =
+      FEEvaluationBaseData<dim, typename Number::value_type, false, Number>;
+
     static void
-    evaluate(const unsigned int                            n_components,
-             const EvaluationFlags::EvaluationFlags        evaluation_flag,
-             const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-             const Number *                                values_dofs,
-             Number *                                      values_quad,
-             Number *                                      gradients_quad,
-             Number *                                      hessians_quad,
-             Number *                                      scratch_data);
+    evaluate(const unsigned int                     n_components,
+             const EvaluationFlags::EvaluationFlags evaluation_flag,
+             const Number *                         values_dofs,
+             EvaluationData &                       eval);
 
     static void
-    integrate(const unsigned int                            n_components,
-              const EvaluationFlags::EvaluationFlags        evaluation_flag,
-              const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-              Number *                                      values_dofs,
-              Number *                                      values_quad,
-              Number *                                      gradients_quad,
-              Number *                                      hessians_quad,
-              Number *                                      scratch_data,
-              const bool add_into_values_array);
+    integrate(const unsigned int                     n_components,
+              const EvaluationFlags::EvaluationFlags evaluation_flag,
+              Number *                               values_dofs,
+              EvaluationData &                       eval,
+              const bool                             add_into_values_array);
   };
 
 
@@ -1623,19 +1644,17 @@ namespace internal
     n_q_points_1d,
     Number>::evaluate(const unsigned int                     n_components,
                       const EvaluationFlags::EvaluationFlags evaluation_flag,
-                      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-                      const Number *                                values_dofs,
-                      Number *                                      values_quad,
-                      Number *gradients_quad,
-                      Number *hessians_quad,
-                      Number *)
+                      const Number *                         values_dofs,
+                      EvaluationData &                       eval)
   {
+    const auto &shape_data = eval.get_shape_info().data[0];
+
     Assert(n_q_points_1d > fe_degree,
            ExcMessage("You lose information when going to a collocation space "
                       "of lower degree, so the evaluation results would be "
                       "wrong. Thus, this class does not permit the desired "
                       "operation."));
-    constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
+    constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
 
     for (unsigned int c = 0; c < n_components; ++c)
       {
@@ -1647,28 +1666,22 @@ namespace internal
           n_q_points_1d,
           Number,
           Number>::do_forward(1,
-                              shape_info.data.front().shape_values_eo,
-                              values_dofs,
-                              values_quad);
+                              shape_data.shape_values_eo,
+                              values_dofs +
+                                c * Utilities::pow(fe_degree + 1, dim),
+                              eval.begin_values() + c * n_q_points);
 
         // apply derivatives in the collocation space
         if (evaluation_flag &
             (EvaluationFlags::gradients | EvaluationFlags::hessians))
-          FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::evaluate(
-            1,
-            evaluation_flag &
-              (EvaluationFlags::gradients | EvaluationFlags::hessians),
-            shape_info,
-            values_quad,
-            nullptr,
-            gradients_quad,
-            hessians_quad,
-            nullptr);
-
-        values_dofs += shape_info.dofs_per_component_on_cell;
-        values_quad += n_q_points;
-        gradients_quad += dim * n_q_points;
-        hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
+          FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::
+            do_evaluate(shape_data,
+                        evaluation_flag & (EvaluationFlags::gradients |
+                                           EvaluationFlags::hessians),
+                        eval.begin_values() + c * n_q_points,
+                        eval.begin_gradients() + c * dim * n_q_points,
+                        eval.begin_hessians() +
+                          c * dim * (dim + 1) / 2 * n_q_points);
       }
   }
 
@@ -1682,24 +1695,18 @@ namespace internal
     n_q_points_1d,
     Number>::integrate(const unsigned int                     n_components,
                        const EvaluationFlags::EvaluationFlags integration_flag,
-                       const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-                       Number *   values_dofs,
-                       Number *   values_quad,
-                       Number *   gradients_quad,
-                       Number *   hessians_quad,
-                       Number *   scratch_data,
+                       Number *                               values_dofs,
+                       EvaluationData &                       eval,
                        const bool add_into_values_array)
   {
+    const auto &shape_data = eval.get_shape_info().data[0];
+
     Assert(n_q_points_1d > fe_degree,
            ExcMessage("You lose information when going to a collocation space "
                       "of lower degree, so the evaluation results would be "
                       "wrong. Thus, this class does not permit the desired "
                       "operation."));
-    AssertDimension(
-      shape_info.data.front().shape_gradients_collocation_eo.size(),
-      (n_q_points_1d + 1) / 2 * n_q_points_1d);
     constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
-    constexpr unsigned int hdim       = (dim * (dim + 1)) / 2;
 
     for (unsigned int c = 0; c < n_components; ++c)
       {
@@ -1707,17 +1714,15 @@ namespace internal
         if (integration_flag &
             (EvaluationFlags::gradients | EvaluationFlags::hessians))
           FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::
-            integrate(1,
-                      integration_flag & (EvaluationFlags::gradients |
-                                          EvaluationFlags::hessians),
-                      shape_info,
-                      values_quad,
-                      nullptr,
-                      gradients_quad,
-                      hessians_quad,
-                      scratch_data,
-                      /*add_into_values_array=*/integration_flag &
-                        EvaluationFlags::values);
+            do_integrate(shape_data,
+                         integration_flag & (EvaluationFlags::gradients |
+                                             EvaluationFlags::hessians),
+                         eval.begin_values() + c * n_q_points,
+                         eval.begin_gradients() + c * dim * n_q_points,
+                         eval.begin_hessians() +
+                           c * dim * (dim + 1) / 2 * n_q_points,
+                         /*add_into_values_array=*/
+                         integration_flag & EvaluationFlags::values);
 
         // transform back to the original space
         FEEvaluationImplBasisChange<
@@ -1728,14 +1733,11 @@ namespace internal
           n_q_points_1d,
           Number,
           Number>::do_backward(1,
-                               shape_info.data.front().shape_values_eo,
+                               shape_data.shape_values_eo,
                                add_into_values_array,
-                               values_quad,
-                               values_dofs);
-        hessians_quad += hdim * n_q_points;
-        gradients_quad += dim * n_q_points;
-        values_quad += n_q_points;
-        values_dofs += shape_info.dofs_per_component_on_cell;
+                               eval.begin_values() + c * n_q_points,
+                               values_dofs +
+                                 c * Utilities::pow(fe_degree + 1, dim));
       }
   }
 
@@ -1761,146 +1763,100 @@ namespace internal
   {
     template <int fe_degree, int n_q_points_1d>
     static bool
-    run(const unsigned int                                      n_components,
-        const EvaluationFlags::EvaluationFlags                  evaluation_flag,
-        const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-        Number *values_dofs_actual,
-        Number *values_quad,
-        Number *gradients_quad,
-        Number *hessians_quad,
-        Number *scratch_data)
+    run(const unsigned int                     n_components,
+        const EvaluationFlags::EvaluationFlags evaluation_flag,
+        const Number *                         values_dofs,
+        FEEvaluationBaseData<dim, typename Number::value_type, false, Number>
+          &eval)
     {
+      const auto element_type = eval.get_shape_info().element_type;
+      using ElementType       = internal::MatrixFreeFunctions::ElementType;
+
       // We enable a transformation to collocation for derivatives if it gives
       // correct results (first condition), if it is the most efficient choice
       // in terms of operation counts (second condition) and if we were able to
       // initialize the fields in shape_info.templates.h from the polynomials
       // (third condition).
-      static constexpr bool use_collocation =
-        n_q_points_1d > fe_degree && n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
-        n_q_points_1d < 200;
-
-      Assert(shape_info.data.size() == 1 ||
-               (shape_info.data.size() == dim &&
-                shape_info.element_type ==
-                  internal::MatrixFreeFunctions::tensor_general),
+      constexpr bool use_collocation = (n_q_points_1d > fe_degree) &&
+                                       (n_q_points_1d < 200) &&
+                                       (n_q_points_1d <= 3 * fe_degree / 2 + 1);
+
+      Assert(eval.get_shape_info().data.size() == 1 ||
+               (eval.get_shape_info().data.size() == dim &&
+                element_type == ElementType::tensor_general),
              ExcNotImplemented());
 
       if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
-          shape_info.element_type ==
-            internal::MatrixFreeFunctions::tensor_symmetric_collocation)
+          element_type == ElementType::tensor_symmetric_collocation)
         {
           internal::FEEvaluationImplCollocation<dim, fe_degree, Number>::
-            evaluate(n_components,
-                     evaluation_flag,
-                     shape_info,
-                     values_dofs_actual,
-                     values_quad,
-                     gradients_quad,
-                     hessians_quad,
-                     scratch_data);
+            evaluate(n_components, evaluation_flag, values_dofs, eval);
         }
       // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
       // shape_info.h for more details
       else if (fe_degree >= 0 && use_collocation &&
-               shape_info.element_type <=
-                 internal::MatrixFreeFunctions::tensor_symmetric)
+               element_type <= ElementType::tensor_symmetric)
         {
           internal::FEEvaluationImplTransformToCollocation<
             dim,
             fe_degree,
             n_q_points_1d,
-            Number>::evaluate(n_components,
-                              evaluation_flag,
-                              shape_info,
-                              values_dofs_actual,
-                              values_quad,
-                              gradients_quad,
-                              hessians_quad,
-                              scratch_data);
+            Number>::evaluate(n_components, evaluation_flag, values_dofs, eval);
         }
-      else if (fe_degree >= 0 &&
-               shape_info.element_type <=
-                 internal::MatrixFreeFunctions::tensor_symmetric)
+      else if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric)
         {
-          internal::FEEvaluationImpl<
-            internal::MatrixFreeFunctions::tensor_symmetric,
-            dim,
-            fe_degree,
-            n_q_points_1d,
-            Number>::evaluate(n_components,
-                              evaluation_flag,
-                              shape_info,
-                              values_dofs_actual,
-                              values_quad,
-                              gradients_quad,
-                              hessians_quad,
-                              scratch_data);
+          internal::FEEvaluationImpl<ElementType::tensor_symmetric,
+                                     dim,
+                                     fe_degree,
+                                     n_q_points_1d,
+                                     Number>::evaluate(n_components,
+                                                       evaluation_flag,
+                                                       values_dofs,
+                                                       eval);
         }
-      else if (shape_info.element_type ==
-               internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+      else if (element_type == ElementType::tensor_symmetric_plus_dg0)
         {
-          internal::FEEvaluationImpl<
-            internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
-            dim,
-            fe_degree,
-            n_q_points_1d,
-            Number>::evaluate(n_components,
-                              evaluation_flag,
-                              shape_info,
-                              values_dofs_actual,
-                              values_quad,
-                              gradients_quad,
-                              hessians_quad,
-                              scratch_data);
+          internal::FEEvaluationImpl<ElementType::tensor_symmetric_plus_dg0,
+                                     dim,
+                                     fe_degree,
+                                     n_q_points_1d,
+                                     Number>::evaluate(n_components,
+                                                       evaluation_flag,
+                                                       values_dofs,
+                                                       eval);
         }
-      else if (shape_info.element_type ==
-               internal::MatrixFreeFunctions::truncated_tensor)
+      else if (element_type == ElementType::truncated_tensor)
         {
-          internal::FEEvaluationImpl<
-            internal::MatrixFreeFunctions::truncated_tensor,
-            dim,
-            fe_degree,
-            n_q_points_1d,
-            Number>::evaluate(n_components,
-                              evaluation_flag,
-                              shape_info,
-                              values_dofs_actual,
-                              values_quad,
-                              gradients_quad,
-                              hessians_quad,
-                              scratch_data);
+          internal::FEEvaluationImpl<ElementType::truncated_tensor,
+                                     dim,
+                                     fe_degree,
+                                     n_q_points_1d,
+                                     Number>::evaluate(n_components,
+                                                       evaluation_flag,
+                                                       values_dofs,
+                                                       eval);
         }
-      else if (shape_info.element_type ==
-               internal::MatrixFreeFunctions::tensor_none)
+      else if (element_type == ElementType::tensor_none)
         {
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_none,
+          internal::FEEvaluationImpl<ElementType::tensor_none,
                                      dim,
                                      fe_degree,
                                      n_q_points_1d,
                                      Number>::evaluate(n_components,
                                                        evaluation_flag,
-                                                       shape_info,
-                                                       values_dofs_actual,
-                                                       values_quad,
-                                                       gradients_quad,
-                                                       hessians_quad,
-                                                       scratch_data);
+                                                       values_dofs,
+                                                       eval);
         }
       else
         {
-          internal::FEEvaluationImpl<
-            internal::MatrixFreeFunctions::tensor_general,
-            dim,
-            fe_degree,
-            n_q_points_1d,
-            Number>::evaluate(n_components,
-                              evaluation_flag,
-                              shape_info,
-                              values_dofs_actual,
-                              values_quad,
-                              gradients_quad,
-                              hessians_quad,
-                              scratch_data);
+          internal::FEEvaluationImpl<ElementType::tensor_general,
+                                     dim,
+                                     fe_degree,
+                                     n_q_points_1d,
+                                     Number>::evaluate(n_components,
+                                                       evaluation_flag,
+                                                       values_dofs,
+                                                       eval);
         }
 
       return false;
@@ -1931,49 +1887,42 @@ namespace internal
     static bool
     run(const unsigned int                     n_components,
         const EvaluationFlags::EvaluationFlags integration_flag,
-        const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-        Number *   values_dofs_actual,
-        Number *   values_quad,
-        Number *   gradients_quad,
-        Number *   hessians_quad,
-        Number *   scratch_data,
+        Number *                               values_dofs,
+        FEEvaluationBaseData<dim, typename Number::value_type, false, Number>
+          &        eval,
         const bool sum_into_values_array)
     {
+      const auto element_type = eval.get_shape_info().element_type;
+      using ElementType       = internal::MatrixFreeFunctions::ElementType;
+
       // We enable a transformation to collocation for derivatives if it gives
       // correct results (first condition), if it is the most efficient choice
       // in terms of operation counts (second condition) and if we were able to
       // initialize the fields in shape_info.templates.h from the polynomials
       // (third condition).
-      constexpr bool use_collocation = n_q_points_1d > fe_degree &&
-                                       n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
-                                       n_q_points_1d < 200;
-
-      Assert(shape_info.data.size() == 1 ||
-               (shape_info.data.size() == dim &&
-                shape_info.element_type ==
-                  internal::MatrixFreeFunctions::tensor_general),
+      constexpr bool use_collocation = (n_q_points_1d > fe_degree) &&
+                                       (n_q_points_1d < 200) &&
+                                       (n_q_points_1d <= 3 * fe_degree / 2 + 1);
+
+      Assert(eval.get_shape_info().data.size() == 1 ||
+               (eval.get_shape_info().data.size() == dim &&
+                element_type == ElementType::tensor_general),
              ExcNotImplemented());
 
       if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
-          shape_info.element_type ==
-            internal::MatrixFreeFunctions::tensor_symmetric_collocation)
+          element_type == ElementType::tensor_symmetric_collocation)
         {
           internal::FEEvaluationImplCollocation<dim, fe_degree, Number>::
             integrate(n_components,
                       integration_flag,
-                      shape_info,
-                      values_dofs_actual,
-                      values_quad,
-                      gradients_quad,
-                      hessians_quad,
-                      scratch_data,
+                      values_dofs,
+                      eval,
                       sum_into_values_array);
         }
       // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
       // shape_info.h for more details
       else if (fe_degree >= 0 && use_collocation &&
-               shape_info.element_type <=
-                 internal::MatrixFreeFunctions::tensor_symmetric)
+               element_type <= ElementType::tensor_symmetric)
         {
           internal::FEEvaluationImplTransformToCollocation<
             dim,
@@ -1981,102 +1930,69 @@ namespace internal
             n_q_points_1d,
             Number>::integrate(n_components,
                                integration_flag,
-                               shape_info,
-                               values_dofs_actual,
-                               values_quad,
-                               gradients_quad,
-                               hessians_quad,
-                               scratch_data,
+                               values_dofs,
+                               eval,
                                sum_into_values_array);
         }
-      else if (fe_degree >= 0 &&
-               shape_info.element_type <=
-                 internal::MatrixFreeFunctions::tensor_symmetric)
+      else if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric)
         {
-          internal::FEEvaluationImpl<
-            internal::MatrixFreeFunctions::tensor_symmetric,
-            dim,
-            fe_degree,
-            n_q_points_1d,
-            Number>::integrate(n_components,
-                               integration_flag,
-                               shape_info,
-                               values_dofs_actual,
-                               values_quad,
-                               gradients_quad,
-                               hessians_quad,
-                               scratch_data,
-                               sum_into_values_array);
+          internal::FEEvaluationImpl<ElementType::tensor_symmetric,
+                                     dim,
+                                     fe_degree,
+                                     n_q_points_1d,
+                                     Number>::integrate(n_components,
+                                                        integration_flag,
+                                                        values_dofs,
+                                                        eval,
+                                                        sum_into_values_array);
         }
-      else if (shape_info.element_type ==
-               internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+      else if (element_type == ElementType::tensor_symmetric_plus_dg0)
         {
-          internal::FEEvaluationImpl<
-            internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
-            dim,
-            fe_degree,
-            n_q_points_1d,
-            Number>::integrate(n_components,
-                               integration_flag,
-                               shape_info,
-                               values_dofs_actual,
-                               values_quad,
-                               gradients_quad,
-                               hessians_quad,
-                               scratch_data,
-                               sum_into_values_array);
+          internal::FEEvaluationImpl<ElementType::tensor_symmetric_plus_dg0,
+                                     dim,
+                                     fe_degree,
+                                     n_q_points_1d,
+                                     Number>::integrate(n_components,
+                                                        integration_flag,
+                                                        values_dofs,
+                                                        eval,
+                                                        sum_into_values_array);
         }
-      else if (shape_info.element_type ==
-               internal::MatrixFreeFunctions::truncated_tensor)
+      else if (element_type == ElementType::truncated_tensor)
         {
-          internal::FEEvaluationImpl<
-            internal::MatrixFreeFunctions::truncated_tensor,
-            dim,
-            fe_degree,
-            n_q_points_1d,
-            Number>::integrate(n_components,
-                               integration_flag,
-                               shape_info,
-                               values_dofs_actual,
-                               values_quad,
-                               gradients_quad,
-                               hessians_quad,
-                               scratch_data,
-                               sum_into_values_array);
+          internal::FEEvaluationImpl<ElementType::truncated_tensor,
+                                     dim,
+                                     fe_degree,
+                                     n_q_points_1d,
+                                     Number>::integrate(n_components,
+                                                        integration_flag,
+                                                        values_dofs,
+                                                        eval,
+                                                        sum_into_values_array);
         }
-      else if (shape_info.element_type ==
-               internal::MatrixFreeFunctions::tensor_none)
+      else if (element_type == ElementType::tensor_none)
         {
-          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_none,
+          internal::FEEvaluationImpl<ElementType::tensor_none,
                                      dim,
                                      fe_degree,
                                      n_q_points_1d,
                                      Number>::integrate(n_components,
                                                         integration_flag,
-                                                        shape_info,
-                                                        values_dofs_actual,
-                                                        values_quad,
-                                                        gradients_quad,
-                                                        hessians_quad,
-                                                        scratch_data,
+                                                        values_dofs,
+                                                        eval,
                                                         sum_into_values_array);
         }
       else
         {
-          internal::FEEvaluationImpl<
-            internal::MatrixFreeFunctions::tensor_general,
-            dim,
-            fe_degree,
-            n_q_points_1d,
-            Number>::integrate(n_components,
-                               integration_flag,
-                               shape_info,
-                               values_dofs_actual,
-                               values_quad,
-                               gradients_quad,
-                               hessians_quad,
-                               scratch_data,
-                               sum_into_values_array);
+          internal::FEEvaluationImpl<ElementType::tensor_general,
+                                     dim,
+                                     fe_degree,
+                                     n_q_points_1d,
+                                     Number>::integrate(n_components,
+                                                        integration_flag,
+                                                        values_dofs,
+                                                        eval,
+                                                        sum_into_values_array);
         }
 
       return false;
@@ -2098,96 +2014,76 @@ namespace internal
     // able to initialize the fields in shape_info.templates.h from the
     // polynomials (fourth condition).
     static constexpr bool use_collocation =
-      symmetric_evaluate &&
-      n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
-      n_q_points_1d < 200;
+      symmetric_evaluate && (n_q_points_1d > fe_degree) &&
+      (n_q_points_1d < 200) && (n_q_points_1d <= 3 * fe_degree / 2 + 1);
+
+    using Eval =
+      EvaluatorTensorProduct<symmetric_evaluate ? internal::evaluate_evenodd :
+                                                  internal::evaluate_general,
+                             dim - 1,
+                             fe_degree + 1,
+                             n_q_points_1d,
+                             Number>;
+
+    using EvaluationData =
+      FEEvaluationBaseData<dim, typename Number::value_type, true, Number>;
+
+    static Eval
+    create_evaluator_tensor_product(
+      const internal::MatrixFreeFunctions::UnivariateShapeData<Number> &data,
+      const unsigned int subface_index,
+      const unsigned int direction)
+    {
+      if (symmetric_evaluate)
+        return Eval(data.shape_values_eo,
+                    data.shape_gradients_eo,
+                    data.shape_hessians_eo,
+                    data.fe_degree + 1,
+                    data.n_q_points_1d);
+      else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
+        return Eval(data.shape_values,
+                    data.shape_gradients,
+                    data.shape_hessians,
+                    data.fe_degree + 1,
+                    data.n_q_points_1d);
+      else
+        {
+          const unsigned int index =
+            direction == 0 ? subface_index % 2 : subface_index / 2;
+          return Eval(data.values_within_subface[index],
+                      data.gradients_within_subface[index],
+                      data.hessians_within_subface[index],
+                      data.fe_degree + 1,
+                      data.n_q_points_1d);
+        }
+    }
 
     static void
-    evaluate_in_face(const unsigned int                            n_components,
-                     const MatrixFreeFunctions::ShapeInfo<Number> &data,
-                     Number *                                      values_dofs,
-                     Number *                                      values_quad,
-                     Number *           gradients_quad,
-                     Number *           hessians_quad,
-                     Number *           scratch_data,
-                     const bool         evaluate_val,
-                     const bool         evaluate_grad,
-                     const bool         evaluate_hessian,
-                     const unsigned int subface_index)
+    evaluate_in_face(
+      const unsigned int                                      n_components,
+      const EvaluationFlags::EvaluationFlags                  evaluation_flag,
+      const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
+      Number *                                                values_dofs,
+      Number *                                                values_quad,
+      Number *                                                gradients_quad,
+      Number *                                                hessians_quad,
+      Number *                                                scratch_data,
+      const unsigned int                                      subface_index)
     {
-      const AlignedVector<Number> &val0 =
-        symmetric_evaluate ?
-          data.data.front().shape_values_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_values :
-             data.data.front().values_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &val1 =
-        symmetric_evaluate ?
-          data.data.front().shape_values_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_values :
-             data.data.front().values_within_subface[subface_index / 2]);
-
-      const AlignedVector<Number> &grad0 =
-        symmetric_evaluate ?
-          data.data.front().shape_gradients_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_gradients :
-             data.data.front().gradients_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &grad1 =
-        symmetric_evaluate ?
-          data.data.front().shape_gradients_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_gradients :
-             data.data.front().gradients_within_subface[subface_index / 2]);
-
-      const AlignedVector<Number> &hessian0 =
-        symmetric_evaluate ?
-          data.data.front().shape_hessians_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_hessians :
-             data.data.front().hessians_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &hessian1 =
-        symmetric_evaluate ?
-          data.data.front().shape_hessians_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_hessians :
-             data.data.front().hessians_within_subface[subface_index / 2]);
-
-      using Eval =
-        internal::EvaluatorTensorProduct<symmetric_evaluate ?
-                                           internal::evaluate_evenodd :
-                                           internal::evaluate_general,
-                                         dim - 1,
-                                         fe_degree + 1,
-                                         n_q_points_1d,
-                                         Number>;
-      Eval eval0(val0,
-                 grad0,
-                 hessian0,
-                 data.data.front().fe_degree + 1,
-                 data.data.front().n_q_points_1d);
-      Eval eval1(val1,
-                 grad1,
-                 hessian1,
-                 data.data.front().fe_degree + 1,
-                 data.data.front().n_q_points_1d);
-
-      const unsigned int size_deg =
-        fe_degree > -1 ?
-          Utilities::pow(fe_degree + 1, dim - 1) :
-          (dim > 1 ?
-             Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
-             1);
+      Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
+      Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
 
-      const unsigned int n_q_points = fe_degree > -1 ?
-                                        Utilities::pow(n_q_points_1d, dim - 1) :
-                                        data.n_q_points_face;
+      const std::size_t size_deg =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         Utilities::pow(data.fe_degree + 1, dim - 1);
+      const std::size_t n_q_points =
+        fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
+                         Utilities::pow(data.n_q_points_1d, dim - 1);
 
       // keep a copy of the original pointer for the case of the Hessians
       Number *values_dofs_ptr = values_dofs;
 
-      if (evaluate_val && !evaluate_grad)
+      if (evaluation_flag == EvaluationFlags::values)
         for (unsigned int c = 0; c < n_components; ++c)
           {
             switch (dim)
@@ -2211,7 +2107,7 @@ namespace internal
             values_dofs += 3 * size_deg;
             values_quad += n_q_points;
           }
-      else if (evaluate_grad)
+      else if (evaluation_flag & EvaluationFlags::gradients)
         for (unsigned int c = 0; c < n_components; ++c)
           {
             switch (dim)
@@ -2229,10 +2125,9 @@ namespace internal
                         n_q_points_1d,
                         n_q_points_1d,
                         Number>
-                        eval_grad(
-                          AlignedVector<Number>(),
-                          data.data.front().shape_gradients_collocation_eo,
-                          AlignedVector<Number>());
+                        eval_grad(AlignedVector<Number>(),
+                                  data.shape_gradients_collocation_eo,
+                                  AlignedVector<Number>());
                       eval_grad.template gradients<0, true, false>(
                         values_quad, gradients_quad);
                       eval_grad.template gradients<1, true, false>(
@@ -2253,7 +2148,7 @@ namespace internal
                                                                gradients_quad +
                                                                  n_q_points);
 
-                      if (evaluate_val == true)
+                      if (evaluation_flag & EvaluationFlags::values)
                         eval1.template values<1, true, false>(scratch_data,
                                                               values_quad);
                     }
@@ -2267,11 +2162,10 @@ namespace internal
                 case 2:
                   eval0.template values<0, true, false>(values_dofs + size_deg,
                                                         gradients_quad +
-                                                          (dim - 1) *
-                                                            n_q_points);
+                                                          n_q_points);
                   eval0.template gradients<0, true, false>(values_dofs,
                                                            gradients_quad);
-                  if (evaluate_val == true)
+                  if (evaluation_flag & EvaluationFlags::values)
                     eval0.template values<0, true, false>(values_dofs,
                                                           values_quad);
                   break;
@@ -2287,10 +2181,9 @@ namespace internal
             gradients_quad += dim * n_q_points;
           }
 
-      if (evaluate_hessian)
+      if (evaluation_flag & EvaluationFlags::hessians)
         {
-          values_dofs                 = values_dofs_ptr;
-          constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
+          values_dofs = values_dofs_ptr;
           for (unsigned int c = 0; c < n_components; ++c)
             {
               switch (dim)
@@ -2359,97 +2252,38 @@ namespace internal
                     AssertThrow(false, ExcNotImplemented());
                 }
               values_dofs += 3 * size_deg;
-              hessians_quad += hdim * n_q_points;
+              hessians_quad += dim * (dim + 1) / 2 * n_q_points;
             }
         }
     }
 
     static void
-    integrate_in_face(const unsigned int n_components,
-                      const MatrixFreeFunctions::ShapeInfo<Number> &data,
-                      Number *                                      values_dofs,
-                      Number *                                      values_quad,
-                      Number *           gradients_quad,
-                      Number *           hessians_quad,
-                      Number *           scratch_data,
-                      const bool         integrate_val,
-                      const bool         integrate_grad,
-                      const bool         integrate_hessian,
-                      const unsigned int subface_index)
+    integrate_in_face(
+      const unsigned int                                      n_components,
+      const EvaluationFlags::EvaluationFlags                  integration_flag,
+      const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
+      Number *                                                values_dofs,
+      Number *                                                values_quad,
+      Number *                                                gradients_quad,
+      Number *                                                hessians_quad,
+      Number *                                                scratch_data,
+      const unsigned int                                      subface_index)
     {
-      const AlignedVector<Number> &val0 =
-        symmetric_evaluate ?
-          data.data.front().shape_values_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_values :
-             data.data.front().values_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &val1 =
-        symmetric_evaluate ?
-          data.data.front().shape_values_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_values :
-             data.data.front().values_within_subface[subface_index / 2]);
-
-      const AlignedVector<Number> &grad0 =
-        symmetric_evaluate ?
-          data.data.front().shape_gradients_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_gradients :
-             data.data.front().gradients_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &grad1 =
-        symmetric_evaluate ?
-          data.data.front().shape_gradients_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_gradients :
-             data.data.front().gradients_within_subface[subface_index / 2]);
-
-      const AlignedVector<Number> &hessian0 =
-        symmetric_evaluate ?
-          data.data.front().shape_hessians_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_hessians :
-             data.data.front().hessians_within_subface[subface_index % 2]);
-      const AlignedVector<Number> &hessian1 =
-        symmetric_evaluate ?
-          data.data.front().shape_hessians_eo :
-          (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.data.front().shape_hessians :
-             data.data.front().hessians_within_subface[subface_index / 2]);
-
-      using Eval =
-        internal::EvaluatorTensorProduct<symmetric_evaluate ?
-                                           internal::evaluate_evenodd :
-                                           internal::evaluate_general,
-                                         dim - 1,
-                                         fe_degree + 1,
-                                         n_q_points_1d,
-                                         Number>;
-      Eval eval0(val0,
-                 grad0,
-                 hessian0,
-                 data.data.front().fe_degree + 1,
-                 data.data.front().n_q_points_1d);
-      Eval eval1(val1,
-                 grad1,
-                 hessian1,
-                 data.data.front().fe_degree + 1,
-                 data.data.front().n_q_points_1d);
-
-      const unsigned int size_deg =
+      Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
+      Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
+
+      const std::size_t size_deg =
         fe_degree > -1 ?
           Utilities::pow(fe_degree + 1, dim - 1) :
-          (dim > 1 ?
-             Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
-             1);
-
-      const unsigned int n_q_points = fe_degree > -1 ?
-                                        Utilities::pow(n_q_points_1d, dim - 1) :
-                                        data.n_q_points_face;
+          (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
+      const std::size_t n_q_points =
+        fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
+                         Utilities::pow(data.n_q_points_1d, dim - 1);
 
       // keep a copy of the original pointer for the case of the Hessians
       Number *values_dofs_ptr = values_dofs;
 
-      if (integrate_val && !integrate_grad)
+      if (integration_flag == EvaluationFlags::values)
         for (unsigned int c = 0; c < n_components; ++c)
           {
             switch (dim)
@@ -2473,7 +2307,7 @@ namespace internal
             values_dofs += 3 * size_deg;
             values_quad += n_q_points;
           }
-      else if (integrate_grad)
+      else if (integration_flag & EvaluationFlags::gradients)
         for (unsigned int c = 0; c < n_components; ++c)
           {
             switch (dim)
@@ -2494,11 +2328,10 @@ namespace internal
                         n_q_points_1d,
                         n_q_points_1d,
                         Number>
-                        eval_grad(
-                          AlignedVector<Number>(),
-                          data.data.front().shape_gradients_collocation_eo,
-                          AlignedVector<Number>());
-                      if (integrate_val)
+                        eval_grad(AlignedVector<Number>(),
+                                  data.shape_gradients_collocation_eo,
+                                  AlignedVector<Number>());
+                      if (integration_flag & EvaluationFlags::values)
                         eval_grad.template gradients<1, false, true>(
                           gradients_quad + n_q_points, values_quad);
                       else
@@ -2513,7 +2346,7 @@ namespace internal
                     }
                   else
                     {
-                      if (integrate_val)
+                      if (integration_flag & EvaluationFlags::values)
                         {
                           eval1.template values<1, false, false>(values_quad,
                                                                  scratch_data);
@@ -2540,7 +2373,7 @@ namespace internal
                     gradients_quad + n_q_points, values_dofs + size_deg);
                   eval0.template gradients<0, false, false>(gradients_quad,
                                                             values_dofs);
-                  if (integrate_val == true)
+                  if (integration_flag & EvaluationFlags::values)
                     eval0.template values<0, false, true>(values_quad,
                                                           values_dofs);
                   break;
@@ -2556,10 +2389,9 @@ namespace internal
             gradients_quad += dim * n_q_points;
           }
 
-      if (integrate_hessian)
+      if (integration_flag & EvaluationFlags::hessians)
         {
-          values_dofs                 = values_dofs_ptr;
-          constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
+          values_dofs = values_dofs_ptr;
           for (unsigned int c = 0; c < n_components; ++c)
             {
               switch (dim)
@@ -2568,7 +2400,8 @@ namespace internal
                     // grad xx
                     eval1.template values<1, false, false>(hessians_quad,
                                                            scratch_data);
-                    if (integrate_val || integrate_grad)
+                    if (integration_flag &
+                        (EvaluationFlags::values | EvaluationFlags::gradients))
                       eval0.template hessians<0, false, true>(scratch_data,
                                                               values_dofs);
                     else
@@ -2601,7 +2434,7 @@ namespace internal
                     eval1.template values<1, false, false>(hessians_quad +
                                                              4 * n_q_points,
                                                            scratch_data);
-                    if (integrate_grad)
+                    if (integration_flag & EvaluationFlags::gradients)
                       eval0.template gradients<0, false, true>(scratch_data,
                                                                values_dofs +
                                                                  size_deg);
@@ -2621,7 +2454,8 @@ namespace internal
                     break;
                   case 2:
                     // grad xx
-                    if (integrate_val || integrate_grad)
+                    if (integration_flag &
+                        (EvaluationFlags::values | EvaluationFlags::gradients))
                       eval0.template hessians<0, false, true>(hessians_quad,
                                                               values_dofs);
                     else
@@ -2632,7 +2466,7 @@ namespace internal
                     eval0.template values<0, false, false>(
                       hessians_quad + n_q_points, values_dofs + 2 * size_deg);
                     // grad xy
-                    if (integrate_grad)
+                    if (integration_flag & EvaluationFlags::gradients)
                       eval0.template gradients<0, false, true>(
                         hessians_quad + 2 * n_q_points, values_dofs + size_deg);
                     else
@@ -2641,16 +2475,16 @@ namespace internal
                     break;
                   case 1:
                     values_dofs[2] = hessians_quad[0];
-                    if (!integrate_val)
+                    if (!(integration_flag & EvaluationFlags::values))
                       values_dofs[0] = 0;
-                    if (!integrate_grad)
+                    if (!(integration_flag & EvaluationFlags::gradients))
                       values_dofs[1] = 0;
                     break;
                   default:
                     AssertThrow(false, ExcNotImplemented());
                 }
               values_dofs += 3 * size_deg;
-              hessians_quad += hdim * n_q_points;
+              hessians_quad += dim * (dim + 1) / 2 * n_q_points;
             }
         }
     }
@@ -2664,15 +2498,13 @@ namespace internal
     template <bool do_evaluate, bool add_into_output>
     static void
     interpolate(const unsigned int                            n_components,
+                const EvaluationFlags::EvaluationFlags        flags,
                 const MatrixFreeFunctions::ShapeInfo<Number> &data,
                 const Number *                                input,
                 Number *                                      output,
-                const bool                                    do_gradients,
-                const bool                                    do_hessians,
                 const unsigned int                            face_no)
     {
-      Assert(static_cast<unsigned int>(fe_degree) ==
-                 data.data.front().fe_degree ||
+      Assert(static_cast<unsigned int>(fe_degree) == data.data[0].fe_degree ||
                fe_degree == -1,
              ExcInternalError());
 
@@ -2680,11 +2512,10 @@ namespace internal
         n_components,
         input,
         output,
-        do_gradients,
-        do_hessians,
+        flags,
         face_no,
-        data.data.front().fe_degree + 1,
-        data.data.front().shape_data_on_face,
+        data.data[0].fe_degree + 1,
+        data.data[0].shape_data_on_face,
         data.dofs_per_component_on_cell,
         3 * data.dofs_per_component_on_face);
     }
@@ -2694,16 +2525,15 @@ namespace internal
      */
     template <bool do_evaluate, bool add_into_output>
     static void
-    interpolate_quadrature(const unsigned int n_components,
+    interpolate_quadrature(const unsigned int                     n_components,
+                           const EvaluationFlags::EvaluationFlags flags,
                            const MatrixFreeFunctions::ShapeInfo<Number> &data,
                            const Number *                                input,
                            Number *                                      output,
-                           const bool         do_gradients,
-                           const bool         do_hessians,
                            const unsigned int face_no)
     {
       Assert(static_cast<unsigned int>(fe_degree + 1) ==
-                 data.data.front().quadrature.size() ||
+                 data.data[0].n_q_points_1d ||
                fe_degree == -1,
              ExcInternalError());
 
@@ -2711,8 +2541,7 @@ namespace internal
         n_components,
         input,
         output,
-        do_gradients,
-        do_hessians,
+        flags,
         face_no,
         data.data.front().quadrature.size(),
         data.data.front().quadrature_data_on_face,
@@ -2723,13 +2552,12 @@ namespace internal
   private:
     template <bool do_evaluate, bool add_into_output, int face_direction = 0>
     static void
-    interpolate_generic(const unsigned int n_components,
-                        const Number *     input,
-                        Number *           output,
-                        const bool         do_gradients,
-                        const bool         do_hessians,
-                        const unsigned int face_no,
-                        const unsigned int n_points_1d,
+    interpolate_generic(const unsigned int                     n_components,
+                        const Number *                         input,
+                        Number *                               output,
+                        const EvaluationFlags::EvaluationFlags flag,
+                        const unsigned int                     face_no,
+                        const unsigned int                     n_points_1d,
                         const std::array<AlignedVector<Number>, 2> &shape_data,
                         const unsigned int dofs_per_component_on_cell,
                         const unsigned int dofs_per_component_on_face)
@@ -2756,13 +2584,13 @@ namespace internal
 
           for (unsigned int c = 0; c < n_components; ++c)
             {
-              if (do_hessians)
+              if (flag & EvaluationFlags::hessians)
                 evalf.template apply_face<face_direction,
                                           do_evaluate,
                                           add_into_output,
                                           2,
                                           lex_faces>(input, output);
-              else if (!do_hessians && do_gradients)
+              else if (flag & EvaluationFlags::gradients)
                 evalf.template apply_face<face_direction,
                                           do_evaluate,
                                           add_into_output,
@@ -2786,8 +2614,7 @@ namespace internal
             n_components,
             input,
             output,
-            do_gradients,
-            do_hessians,
+            flag,
             face_no,
             n_points_1d,
             shape_data,
@@ -2906,25 +2733,20 @@ namespace internal
 
   template <typename Number>
   void
-  adjust_for_face_orientation(const unsigned int            dim,
-                              const unsigned int            n_components,
-                              const unsigned int            face_orientation,
-                              const Table<2, unsigned int> &orientation_map,
-                              const bool                    integrate,
-                              const bool                    values,
-                              const bool                    gradients,
-                              const bool                    hessians,
-                              const unsigned int            n_q_points,
-                              Number *                      tmp_values,
-                              Number *                      values_quad,
-                              Number *                      gradients_quad,
-                              Number *                      hessians_quad)
+  adjust_for_face_orientation(const unsigned int dim,
+                              const unsigned int n_components,
+                              const EvaluationFlags::EvaluationFlags flag,
+                              const unsigned int *orientation,
+                              const bool          integrate,
+                              const std::size_t   n_q_points,
+                              Number *            tmp_values,
+                              Number *            values_quad,
+                              Number *            gradients_quad,
+                              Number *            hessians_quad)
   {
-    Assert(face_orientation, ExcInternalError());
-    const unsigned int *orientation = &orientation_map[face_orientation][0];
     for (unsigned int c = 0; c < n_components; ++c)
       {
-        if (values == true)
+        if (flag & EvaluationFlags::values)
           {
             if (integrate)
               for (unsigned int q = 0; q < n_q_points; ++q)
@@ -2935,7 +2757,7 @@ namespace internal
             for (unsigned int q = 0; q < n_q_points; ++q)
               values_quad[c * n_q_points + q] = tmp_values[q];
           }
-        if (gradients == true)
+        if (flag & EvaluationFlags::gradients)
           for (unsigned int d = 0; d < dim; ++d)
             {
               if (integrate)
@@ -2949,7 +2771,7 @@ namespace internal
               for (unsigned int q = 0; q < n_q_points; ++q)
                 gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
             }
-        if (hessians == true)
+        if (flag & EvaluationFlags::hessians)
           {
             const unsigned int hdim = (dim * (dim + 1)) / 2;
             for (unsigned int d = 0; d < hdim; ++d)
@@ -2972,46 +2794,37 @@ namespace internal
 
 
 
-  template <int dim, typename VectorizedArrayType>
+  template <int dim, typename Number>
   struct FEFaceEvaluationImplEvaluateSelector
   {
     template <int fe_degree, int n_q_points_1d>
     static bool
-    run(const unsigned int                                         n_components,
-        const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-        const VectorizedArrayType *                                values_array,
-        VectorizedArrayType *                                      values_quad,
-        VectorizedArrayType *         gradients_quad,
-        VectorizedArrayType *         hessians_quad,
-        VectorizedArrayType *         scratch_data,
-        const bool                    evaluate_values,
-        const bool                    evaluate_gradients,
-        const bool                    evaluate_hessians,
-        const unsigned int            face_no,
-        const unsigned int            subface_index,
-        const unsigned int            face_orientation,
-        const Table<2, unsigned int> &orientation_map)
+    run(const unsigned int                     n_components,
+        const EvaluationFlags::EvaluationFlags evaluation_flag,
+        const Number *                         values_dofs,
+        FEEvaluationBaseData<dim, typename Number::value_type, true, Number>
+          &eval)
     {
+      const auto &       data             = eval.get_shape_info();
+      const auto &       shape_data       = data.data.front();
+      const unsigned int face_no          = eval.get_face_no();
+      const unsigned int face_orientation = eval.get_face_orientation();
+
       if (data.element_type == MatrixFreeFunctions::tensor_none)
         {
           const unsigned int n_dofs     = data.dofs_per_component_on_cell;
           const unsigned int n_q_points = data.n_q_points_faces[face_no];
-          const auto         shape_info = data.data.front();
 
-          using Eval = EvaluatorTensorProduct<evaluate_general,
-                                              1,
-                                              0,
-                                              0,
-                                              VectorizedArrayType,
-                                              VectorizedArrayType>;
+          using Eval =
+            EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
 
-          if (evaluate_values)
+          if (evaluation_flag & EvaluationFlags::values)
             {
               const auto shape_values =
-                &shape_info.shape_values_face(face_no, face_orientation, 0);
+                &shape_data.shape_values_face(face_no, face_orientation, 0);
 
-              auto values_quad_ptr        = values_quad;
-              auto values_dofs_actual_ptr = values_array;
+              auto values_quad_ptr        = eval.begin_values();
+              auto values_dofs_actual_ptr = values_dofs;
 
               Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
               for (unsigned int c = 0; c < n_components; ++c)
@@ -3024,14 +2837,14 @@ namespace internal
                 }
             }
 
-          if (evaluate_gradients)
+          if (evaluation_flag & EvaluationFlags::gradients)
             {
-              auto gradients_quad_ptr     = gradients_quad;
-              auto values_dofs_actual_ptr = values_array;
+              auto gradients_quad_ptr     = eval.begin_gradients();
+              auto values_dofs_actual_ptr = values_dofs;
 
-              std::array<const VectorizedArrayType *, dim> shape_gradients;
+              std::array<const Number *, dim> shape_gradients;
               for (unsigned int d = 0; d < dim; ++d)
-                shape_gradients[d] = &shape_info.shape_gradients_face(
+                shape_gradients[d] = &shape_data.shape_gradients_face(
                   face_no, face_orientation, d, 0);
 
               for (unsigned int c = 0; c < n_components; ++c)
@@ -3053,90 +2866,70 @@ namespace internal
                 }
             }
 
-          if (evaluate_hessians)
-            {
-              AssertThrow(false, ExcNotImplemented());
-            }
+          Assert(!(evaluation_flag & EvaluationFlags::hessians),
+                 ExcNotImplemented());
 
           return true;
         }
 
-      constexpr unsigned int static_dofs_per_face =
-        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                         numbers::invalid_unsigned_int;
       const unsigned int dofs_per_face =
-        fe_degree > -1 ?
-          static_dofs_per_face :
-          Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         Utilities::pow(shape_data.fe_degree + 1, dim - 1);
 
-      VectorizedArrayType *temp1 = scratch_data;
+      Number *temp         = eval.get_scratch_data().begin();
+      Number *scratch_data = temp + 3 * n_components * dofs_per_face;
 
-      FEFaceNormalEvaluationImpl<dim, fe_degree, VectorizedArrayType>::
-        template interpolate<true, false>(n_components,
-                                          data,
-                                          values_array,
-                                          temp1,
-                                          evaluate_gradients,
-                                          evaluate_hessians,
-                                          face_no);
+      FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::template interpolate<
+        true,
+        false>(n_components, evaluation_flag, data, values_dofs, temp, face_no);
 
-      const unsigned int n_q_points_1d_actual =
+      const unsigned int     subface_index = eval.get_subface_index();
+      constexpr unsigned int n_q_points_1d_actual =
         fe_degree > -1 ? n_q_points_1d : 0;
+
       if (fe_degree > -1 &&
           subface_index >= GeometryInfo<dim>::max_children_per_cell &&
           data.element_type <= MatrixFreeFunctions::tensor_symmetric)
-        FEFaceEvaluationImpl<
-          true,
-          dim,
-          fe_degree,
-          n_q_points_1d_actual,
-          VectorizedArrayType>::evaluate_in_face(n_components,
-                                                 data,
-                                                 temp1,
-                                                 values_quad,
-                                                 gradients_quad,
-                                                 hessians_quad,
-                                                 scratch_data + 3 *
-                                                                  n_components *
-                                                                  dofs_per_face,
-                                                 evaluate_values,
-                                                 evaluate_gradients,
-                                                 evaluate_hessians,
-                                                 subface_index);
+        FEFaceEvaluationImpl<true,
+                             dim,
+                             fe_degree,
+                             n_q_points_1d_actual,
+                             Number>::evaluate_in_face(n_components,
+                                                       evaluation_flag,
+                                                       shape_data,
+                                                       temp,
+                                                       eval.begin_values(),
+                                                       eval.begin_gradients(),
+                                                       eval.begin_hessians(),
+                                                       scratch_data,
+                                                       subface_index);
       else
-        FEFaceEvaluationImpl<
-          false,
-          dim,
-          fe_degree,
-          n_q_points_1d_actual,
-          VectorizedArrayType>::evaluate_in_face(n_components,
-                                                 data,
-                                                 temp1,
-                                                 values_quad,
-                                                 gradients_quad,
-                                                 hessians_quad,
-                                                 scratch_data + 3 *
-                                                                  n_components *
-                                                                  dofs_per_face,
-                                                 evaluate_values,
-                                                 evaluate_gradients,
-                                                 evaluate_hessians,
-                                                 subface_index);
+        FEFaceEvaluationImpl<false,
+                             dim,
+                             fe_degree,
+                             n_q_points_1d_actual,
+                             Number>::evaluate_in_face(n_components,
+                                                       evaluation_flag,
+                                                       shape_data,
+                                                       temp,
+                                                       eval.begin_values(),
+                                                       eval.begin_gradients(),
+                                                       eval.begin_hessians(),
+                                                       scratch_data,
+                                                       subface_index);
 
       if (face_orientation)
-        adjust_for_face_orientation(dim,
-                                    n_components,
-                                    face_orientation,
-                                    orientation_map,
-                                    false,
-                                    evaluate_values,
-                                    evaluate_gradients,
-                                    evaluate_hessians,
-                                    data.n_q_points_face,
-                                    scratch_data,
-                                    values_quad,
-                                    gradients_quad,
-                                    hessians_quad);
+        adjust_for_face_orientation(
+          dim,
+          n_components,
+          evaluation_flag,
+          &eval.get_orientation_map()(face_orientation, 0),
+          false,
+          data.n_q_points_face,
+          temp,
+          eval.begin_values(),
+          eval.begin_gradients(),
+          eval.begin_hessians());
 
       return false;
     }
@@ -3144,46 +2937,37 @@ namespace internal
 
 
 
-  template <int dim, typename VectorizedArrayType>
+  template <int dim, typename Number>
   struct FEFaceEvaluationImplIntegrateSelector
   {
     template <int fe_degree, int n_q_points_1d>
     static bool
-    run(const unsigned int                                         n_components,
-        const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-        VectorizedArrayType *                                      values_array,
-        VectorizedArrayType *                                      values_quad,
-        VectorizedArrayType *         gradients_quad,
-        VectorizedArrayType *         hessians_quad,
-        VectorizedArrayType *         scratch_data,
-        const bool                    integrate_values,
-        const bool                    integrate_gradients,
-        const bool                    integrate_hessians,
-        const unsigned int            face_no,
-        const unsigned int            subface_index,
-        const unsigned int            face_orientation,
-        const Table<2, unsigned int> &orientation_map)
+    run(const unsigned int                     n_components,
+        const EvaluationFlags::EvaluationFlags integration_flag,
+        Number *                               values_dofs,
+        FEEvaluationBaseData<dim, typename Number::value_type, true, Number>
+          &eval)
     {
+      const auto &       data             = eval.get_shape_info();
+      const auto &       shape_data       = data.data.front();
+      const unsigned int face_no          = eval.get_face_no();
+      const unsigned int face_orientation = eval.get_face_orientation();
+
       if (data.element_type == MatrixFreeFunctions::tensor_none)
         {
           const unsigned int n_dofs     = data.dofs_per_component_on_cell;
           const unsigned int n_q_points = data.n_q_points_faces[face_no];
-          const auto         shape_info = data.data.front();
 
-          using Eval = EvaluatorTensorProduct<evaluate_general,
-                                              1,
-                                              0,
-                                              0,
-                                              VectorizedArrayType,
-                                              VectorizedArrayType>;
+          using Eval =
+            EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
 
-          if (integrate_values)
+          if (integration_flag & EvaluationFlags::values)
             {
               const auto shape_values =
-                &shape_info.shape_values_face(face_no, face_orientation, 0);
+                &shape_data.shape_values_face(face_no, face_orientation, 0);
 
-              auto values_quad_ptr        = values_quad;
-              auto values_dofs_actual_ptr = values_array;
+              auto values_quad_ptr        = eval.begin_values();
+              auto values_dofs_actual_ptr = values_dofs;
 
               Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
               for (unsigned int c = 0; c < n_components; ++c)
@@ -3196,14 +2980,14 @@ namespace internal
                 }
             }
 
-          if (integrate_gradients)
+          if (integration_flag & EvaluationFlags::gradients)
             {
-              auto gradients_quad_ptr     = gradients_quad;
-              auto values_dofs_actual_ptr = values_array;
+              auto gradients_quad_ptr     = eval.begin_gradients();
+              auto values_dofs_actual_ptr = values_dofs;
 
-              std::array<const VectorizedArrayType *, dim> shape_gradients;
+              std::array<const Number *, dim> shape_gradients;
               for (unsigned int d = 0; d < dim; ++d)
-                shape_gradients[d] = &shape_info.shape_gradients_face(
+                shape_gradients[d] = &shape_data.shape_gradients_face(
                   face_no, face_orientation, d, 0);
 
               for (unsigned int c = 0; c < n_components; ++c)
@@ -3216,7 +3000,8 @@ namespace internal
                                 n_dofs,
                                 n_q_points);
 
-                      if ((integrate_values == false) && d == 0)
+                      if (!(integration_flag & EvaluationFlags::values) &&
+                          d == 0)
                         eval.template gradients<0, false, false>(
                           gradients_quad_ptr, values_dofs_actual_ptr);
                       else
@@ -3229,209 +3014,130 @@ namespace internal
                 }
             }
 
-          if (integrate_hessians)
-            {
-              AssertThrow(false, ExcNotImplemented());
-            }
+          Assert(!(integration_flag & EvaluationFlags::hessians),
+                 ExcNotImplemented());
 
           return true;
         }
 
-      if (face_orientation)
-        adjust_for_face_orientation(dim,
-                                    n_components,
-                                    face_orientation,
-                                    orientation_map,
-                                    true,
-                                    integrate_values,
-                                    integrate_gradients,
-                                    integrate_hessians,
-                                    data.n_q_points_face,
-                                    scratch_data,
-                                    values_quad,
-                                    gradients_quad,
-                                    hessians_quad);
-
-      constexpr unsigned int static_dofs_per_face =
-        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                         numbers::invalid_unsigned_int;
       const unsigned int dofs_per_face =
-        fe_degree > -1 ?
-          static_dofs_per_face :
-          Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         Utilities::pow(shape_data.fe_degree + 1, dim - 1);
+
+      Number *temp         = eval.get_scratch_data().begin();
+      Number *scratch_data = temp + 3 * n_components * dofs_per_face;
 
-      VectorizedArrayType *temp1 = scratch_data;
+      if (face_orientation)
+        adjust_for_face_orientation(
+          dim,
+          n_components,
+          integration_flag,
+          &eval.get_orientation_map()(face_orientation, 0),
+          true,
+          data.n_q_points_face,
+          temp,
+          eval.begin_values(),
+          eval.begin_gradients(),
+          eval.begin_hessians());
 
       const unsigned int n_q_points_1d_actual =
         fe_degree > -1 ? n_q_points_1d : 0;
+      const unsigned int subface_index = eval.get_subface_index();
+
       if (fe_degree > -1 &&
-          subface_index >= GeometryInfo<dim - 1>::max_children_per_cell &&
+          eval.get_subface_index() >=
+            GeometryInfo<dim - 1>::max_children_per_cell &&
           data.element_type <= MatrixFreeFunctions::tensor_symmetric)
-        FEFaceEvaluationImpl<
-          true,
-          dim,
-          fe_degree,
-          n_q_points_1d_actual,
-          VectorizedArrayType>::integrate_in_face(n_components,
-                                                  data,
-                                                  temp1,
-                                                  values_quad,
-                                                  gradients_quad,
-                                                  hessians_quad,
-                                                  scratch_data +
-                                                    3 * n_components *
-                                                      dofs_per_face,
-                                                  integrate_values,
-                                                  integrate_gradients,
-                                                  integrate_hessians,
-                                                  subface_index);
+        FEFaceEvaluationImpl<true,
+                             dim,
+                             fe_degree,
+                             n_q_points_1d_actual,
+                             Number>::integrate_in_face(n_components,
+                                                        integration_flag,
+                                                        shape_data,
+                                                        temp,
+                                                        eval.begin_values(),
+                                                        eval.begin_gradients(),
+                                                        eval.begin_hessians(),
+                                                        scratch_data,
+                                                        subface_index);
       else
-        FEFaceEvaluationImpl<
-          false,
-          dim,
-          fe_degree,
-          n_q_points_1d_actual,
-          VectorizedArrayType>::integrate_in_face(n_components,
-                                                  data,
-                                                  temp1,
-                                                  values_quad,
-                                                  gradients_quad,
-                                                  hessians_quad,
-                                                  scratch_data +
-                                                    3 * n_components *
-                                                      dofs_per_face,
-                                                  integrate_values,
-                                                  integrate_gradients,
-                                                  integrate_hessians,
-                                                  subface_index);
+        FEFaceEvaluationImpl<false,
+                             dim,
+                             fe_degree,
+                             n_q_points_1d_actual,
+                             Number>::integrate_in_face(n_components,
+                                                        integration_flag,
+                                                        shape_data,
+                                                        temp,
+                                                        eval.begin_values(),
+                                                        eval.begin_gradients(),
+                                                        eval.begin_hessians(),
+                                                        scratch_data,
+                                                        subface_index);
 
-      FEFaceNormalEvaluationImpl<dim, fe_degree, VectorizedArrayType>::
-        template interpolate<false, false>(n_components,
-                                           data,
-                                           temp1,
-                                           values_array,
-                                           integrate_gradients,
-                                           integrate_hessians,
-                                           face_no);
+      FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
+        template interpolate<false, false>(
+          n_components, integration_flag, data, temp, values_dofs, face_no);
       return false;
     }
   };
 
 
 
-  template <int n_face_orientations, typename Processor>
-  static bool
-  fe_face_evaluation_process_and_io(Processor &proc)
+  template <int n_face_orientations,
+            typename Processor,
+            typename EvaluationData>
+  void
+  fe_face_evaluation_process_and_io(
+    const Processor &                      proc,
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags evaluation_flag,
+    typename Processor::Number2_ *         src_ptr,
+    const std::vector<ArrayView<const typename Processor::Number_>> *sm_ptr,
+    const EvaluationData &                                           eval,
+    typename Processor::VectorizedArrayType_ *                       temp1)
   {
-    auto  n_components             = proc.n_components;
-    auto  integrate                = proc.integrate;
-    auto  global_vector_ptr        = proc.global_vector_ptr;
-    auto &sm_ptr                   = proc.sm_ptr;
-    auto &data                     = proc.data;
-    auto &dof_info                 = proc.dof_info;
-    auto  values_quad              = proc.values_quad;
-    auto  gradients_quad           = proc.gradients_quad;
-    auto  hessians_quad            = proc.hessians_quad;
-    auto  scratch_data             = proc.scratch_data;
-    auto  do_values                = proc.do_values;
-    auto  do_gradients             = proc.do_gradients;
-    auto  do_hessians              = proc.do_hessians;
-    auto  active_fe_index          = proc.active_fe_index;
-    auto  first_selected_component = proc.first_selected_component;
-    auto  cells                    = proc.cells;
-    auto  face_nos                 = proc.face_nos;
-    auto  subface_index            = proc.subface_index;
-    auto  dof_access_index         = proc.dof_access_index;
-    auto  face_orientations        = proc.face_orientations;
-    auto &orientation_map          = proc.orientation_map;
-
-    static const int dim       = Processor::dim_;
-    static const int fe_degree = Processor::fe_degree_;
-    using VectorizedArrayType  = typename Processor::VectorizedArrayType_;
+    (void)proc;
+    (void)n_components;
+    (void)evaluation_flag;
+    (void)src_ptr;
+    (void)sm_ptr;
+    (void)eval;
+    (void)temp1;
+
+    /*
+    constexpr int dim         = Processor::dim_;
+    constexpr int fe_degree   = Processor::fe_degree_;
+    using VectorizedArrayType = typename Processor::VectorizedArrayType_;
 
     using Number   = typename Processor::Number_;
     using Number2_ = typename Processor::Number2_;
 
-    const unsigned int cell = cells[0];
-
-    // In the case of integration, we do not need to reshuffle the
-    // data at the quadrature points to adjust for the face
-    // orientation if the shape functions are nodal at the cell
-    // boundaries (and we only requested the integration of the
-    // values) or Hermite shape functions are used. These cases are
-    // handled later when the values are written back into the
-    // glrobal vector.
-    if (integrate &&
-        (face_orientations[0] > 0 &&
-         (subface_index < GeometryInfo<dim>::max_children_per_cell ||
-          !(((do_gradients == false &&
-              data.data.front().nodal_at_cell_boundaries == true &&
-              fe_degree > 0) ||
-             (data.element_type ==
-                MatrixFreeFunctions::tensor_symmetric_hermite &&
-              fe_degree > 1)) &&
-            (dof_info.index_storage_variants[dof_access_index][cell] ==
-               MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                 interleaved_contiguous ||
-             dof_info.index_storage_variants[dof_access_index][cell] ==
-               MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                 interleaved_contiguous_strided ||
-             dof_info.index_storage_variants[dof_access_index][cell] ==
-               MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                 interleaved_contiguous_mixed_strides ||
-             dof_info.index_storage_variants[dof_access_index][cell] ==
-               MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                 contiguous)))))
-      {
-        AssertDimension(n_face_orientations, 1);
-        adjust_for_face_orientation(dim,
-                                    n_components,
-                                    face_orientations[0],
-                                    orientation_map,
-                                    true,
-                                    do_values,
-                                    do_gradients,
-                                    do_hessians,
-                                    data.n_q_points_face,
-                                    scratch_data,
-                                    values_quad,
-                                    gradients_quad,
-                                    hessians_quad);
-      }
+    const auto &       shape_data   = eval.get_shape_info().data[0];
+    const auto         element_type = eval.get_shape_info().element_type;
+    constexpr bool     integrate    = Processor::do_integrate;
+    const unsigned int face_no      = eval.get_face_no();
 
     // 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.
     VectorizedArrayType grad_weight =
-      (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
-       data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
-        data.data.front()
-          .shape_data_on_face[0][fe_degree + (integrate ?
-                                                (2 - (face_nos[0] % 2)) :
-                                                (1 + (face_nos[0] % 2)))] :
-        VectorizedArrayType(0.0 /*dummy*/);
+      shape_data
+        .shape_data_on_face[0][fe_degree + integrate ? (2 - face_no % 2) :
+                                                       (1 + face_no % 2)];
 
     constexpr unsigned int static_dofs_per_component =
-      fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
-                       numbers::invalid_unsigned_int;
-    constexpr unsigned int static_dofs_per_face =
-      fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                       numbers::invalid_unsigned_int;
-    const unsigned int dofs_per_face =
-      fe_degree > -1 ? static_dofs_per_face :
-                       Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
-
-    VectorizedArrayType *temp1 = scratch_data;
-
-    const unsigned int dummy = 0;
+      Utilities::pow(fe_degree + 1, dim);
+    constexpr unsigned int dofs_per_face =
+      Utilities::pow(fe_degree + 1, dim - 1);
 
     // re-orientation
     std::array<const unsigned int *, n_face_orientations> orientation = {};
 
     if (n_face_orientations == 1)
-      orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ?
-                         &data.face_orientations[face_orientations[0]][0] :
-                         &dummy;
+      orientation[0] =
+        &eval.get_orientation_map()(eval.get_face_orientation(), 0);
     else
       {
         for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
@@ -3443,7 +3149,7 @@ namespace internal
               continue;
 
             orientation[v] =
-              (data.data.front().nodal_at_cell_boundaries == true) ?
+              (shape_data.nodal_at_cell_boundaries == true) ?
                 &data.face_orientations[face_orientations[v]][0] :
                 &dummy;
           }
@@ -3517,541 +3223,471 @@ namespace internal
 
     for (unsigned int comp = 0; comp < n_components; ++comp)
       {
-        if (integrate)
-          proc.in_face_operation(temp1, comp);
-
         // we can only use the fast functions if we know the polynomial degree
         // as a template parameter (fe_degree != -1), and it only makes sense
         // to use the functions for at least linear functions for values on
         // the faces and quadratic functions for gradients on the faces, so
         // include the switch here
-        if (((do_gradients == false &&
-              data.data.front().nodal_at_cell_boundaries == true &&
-              fe_degree > 0) ||
-             (data.element_type ==
-                MatrixFreeFunctions::tensor_symmetric_hermite &&
-              fe_degree > 1)) &&
-            do_hessians == false)
+        // case 1: contiguous and interleaved indices
+        if (n_face_orientations == 1 &&
+            dof_info.index_storage_variants[dof_access_index][cell] ==
+              MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                interleaved_contiguous)
           {
-            // case 1: contiguous and interleaved indices
-            if (n_face_orientations == 1 &&
-                dof_info.index_storage_variants[dof_access_index][cell] ==
-                  MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                    interleaved_contiguous)
+            AssertDimension(n_face_orientations, 1);
+
+            AssertDimension(
+              dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+              VectorizedArrayType::size());
+            Number2_ *vector_ptr =
+              global_vector_ptr +
+              dof_info
+                .dof_indices_contiguous[dof_access_index]
+                                       [cell * VectorizedArrayType::size()] +
+              (dof_info.component_dof_indices_offset[active_fe_index]
+                                                    [first_selected_component] +
+               comp * static_dofs_per_component) *
+                VectorizedArrayType::size();
+
+            if (fe_degree > 1 && do_gradients == true)
               {
-                AssertDimension(n_face_orientations, 1);
-
-                AssertDimension(
-                  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
-                  VectorizedArrayType::size());
-                Number2_ *vector_ptr =
-                  global_vector_ptr +
-                  dof_info.dof_indices_contiguous[dof_access_index]
-                                                 [cell *
-                                                  VectorizedArrayType::size()] +
-                  (dof_info
-                     .component_dof_indices_offset[active_fe_index]
-                                                  [first_selected_component] +
-                   comp * static_dofs_per_component) *
-                    VectorizedArrayType::size();
-
-                if (fe_degree > 1 && do_gradients == true)
+                for (unsigned int i = 0; i < dofs_per_face; ++i)
                   {
-                    for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    if (n_face_orientations == 1)
                       {
-                        if (n_face_orientations == 1)
-                          {
-                            const unsigned int ind1 =
-                              index_array_hermite[0][2 * i];
-                            const unsigned int ind2 =
-                              index_array_hermite[0][2 * i + 1];
-                            AssertIndexRange(ind1,
-                                             data.dofs_per_component_on_cell);
-                            AssertIndexRange(ind2,
-                                             data.dofs_per_component_on_cell);
-                            const unsigned int i_ = reorientate(0, i);
-                            proc.hermite_grad_vectorized(
-                              temp1[i_],
-                              temp1[i_ + dofs_per_face],
-                              vector_ptr + ind1 * VectorizedArrayType::size(),
-                              vector_ptr + ind2 * VectorizedArrayType::size(),
-                              grad_weight);
-                          }
-                        else
-                          {
-                            Assert(false, ExcNotImplemented());
-                          }
+                        const unsigned int ind1 = index_array_hermite[0][2 * i];
+                        const unsigned int ind2 =
+                          index_array_hermite[0][2 * i + 1];
+                        AssertIndexRange(ind1, data.dofs_per_component_on_cell);
+                        AssertIndexRange(ind2, data.dofs_per_component_on_cell);
+                        const unsigned int i_ = reorientate(0, i);
+                        proc.hermite_grad_vectorized(
+                          temp1[i_],
+                          temp1[i_ + dofs_per_face],
+                          vector_ptr + ind1 * VectorizedArrayType::size(),
+                          vector_ptr + ind2 * VectorizedArrayType::size(),
+                          grad_weight);
                       }
-                  }
-                else
-                  {
-                    for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    else
                       {
-                        if (n_face_orientations == 1)
-                          {
-                            const unsigned int i_  = reorientate(0, i);
-                            const unsigned int ind = index_array_nodal[0][i];
-                            proc.value_vectorized(
-                              temp1[i_],
-                              vector_ptr + ind * VectorizedArrayType::size());
-                          }
-                        else
-                          {
-                            Assert(false, ExcNotImplemented());
-                          }
+                        Assert(false, ExcNotImplemented());
                       }
                   }
               }
-
-            // case 2: contiguous and interleaved indices with fixed stride
-            else if (n_face_orientations == 1 &&
-                     dof_info.index_storage_variants[dof_access_index][cell] ==
-                       MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                         interleaved_contiguous_strided)
+            else
               {
-                AssertDimension(n_face_orientations, 1);
-
-                AssertDimension(
-                  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
-                  VectorizedArrayType::size());
-                const unsigned int *indices =
-                  &dof_info.dof_indices_contiguous[dof_access_index]
-                                                  [cell *
-                                                   VectorizedArrayType::size()];
-                Number2_ *vector_ptr =
-                  global_vector_ptr +
-                  (comp * static_dofs_per_component +
-                   dof_info
-                     .component_dof_indices_offset[active_fe_index]
-                                                  [first_selected_component]) *
-                    VectorizedArrayType::size();
-                if (fe_degree > 1 && do_gradients == true)
+                for (unsigned int i = 0; i < dofs_per_face; ++i)
                   {
-                    for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    if (n_face_orientations == 1)
                       {
-                        if (n_face_orientations == 1)
-                          {
-                            const unsigned int i_ = reorientate(0, i);
-                            const unsigned int ind1 =
-                              index_array_hermite[0][2 * i] *
-                              VectorizedArrayType::size();
-                            const unsigned int ind2 =
-                              index_array_hermite[0][2 * i + 1] *
-                              VectorizedArrayType::size();
-                            proc.hermite_grad_vectorized_indexed(
-                              temp1[i_],
-                              temp1[i_ + dofs_per_face],
-                              vector_ptr + ind1,
-                              vector_ptr + ind2,
-                              grad_weight,
-                              indices,
-                              indices);
-                          }
-                        else
-                          {
-                            Assert(false, ExcNotImplemented());
-                          }
+                        const unsigned int i_  = reorientate(0, i);
+                        const unsigned int ind = index_array_nodal[0][i];
+                        proc.value_vectorized(temp1[i_],
+                                              vector_ptr +
+                                                ind *
+                                                  VectorizedArrayType::size());
                       }
-                  }
-                else
-                  {
-                    for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    else
                       {
-                        if (n_face_orientations == 1)
-                          {
-                            const unsigned int i_ = reorientate(0, i);
-                            const unsigned int ind =
-                              index_array_nodal[0][i] *
-                              VectorizedArrayType::size();
-                            proc.value_vectorized_indexed(temp1[i_],
-                                                          vector_ptr + ind,
-                                                          indices);
-                          }
-                        else
-                          {
-                            Assert(false, ExcNotImplemented());
-                          }
+                        Assert(false, ExcNotImplemented());
                       }
                   }
               }
+          }
 
-            // case 3: contiguous and interleaved indices with mixed stride
-            else if (n_face_orientations == 1 &&
-                     dof_info.index_storage_variants[dof_access_index][cell] ==
-                       MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                         interleaved_contiguous_mixed_strides)
+        // case 2: contiguous and interleaved indices with fixed stride
+        else if (n_face_orientations == 1 &&
+                 dof_info.index_storage_variants[dof_access_index][cell] ==
+                   MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                     interleaved_contiguous_strided)
+          {
+            AssertDimension(n_face_orientations, 1);
+
+            AssertDimension(
+              dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+              VectorizedArrayType::size());
+            const unsigned int *indices =
+              &dof_info
+                 .dof_indices_contiguous[dof_access_index]
+                                        [cell * VectorizedArrayType::size()];
+            Number2_ *vector_ptr =
+              global_vector_ptr +
+              (comp * static_dofs_per_component +
+               dof_info
+                 .component_dof_indices_offset[active_fe_index]
+                                              [first_selected_component]) *
+                VectorizedArrayType::size();
+            if (fe_degree > 1 && do_gradients == true)
               {
-                AssertDimension(n_face_orientations, 1);
-
-                const unsigned int *strides =
-                  &dof_info.dof_indices_interleave_strides
-                     [dof_access_index][cell * VectorizedArrayType::size()];
-                unsigned int indices[VectorizedArrayType::size()];
-                for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
-                  indices[v] =
-                    dof_info.dof_indices_contiguous
-                      [dof_access_index]
-                      [cell * VectorizedArrayType::size() + v] +
-                    (dof_info
-                       .component_dof_indices_offset[active_fe_index]
-                                                    [first_selected_component] +
-                     comp * static_dofs_per_component) *
-                      strides[v];
-                const unsigned int n_filled_lanes =
-                  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
-
-                if (fe_degree > 1 && do_gradients == true)
+                for (unsigned int i = 0; i < dofs_per_face; ++i)
                   {
-                    if (n_filled_lanes == VectorizedArrayType::size())
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
-                        {
-                          if (n_face_orientations == 1)
-                            {
-                              const unsigned int i_ = reorientate(0, i);
-                              unsigned int ind1[VectorizedArrayType::size()];
-                              DEAL_II_OPENMP_SIMD_PRAGMA
-                              for (unsigned int v = 0;
-                                   v < VectorizedArrayType::size();
-                                   ++v)
-                                ind1[v] =
-                                  indices[v] +
-                                  index_array_hermite[0 /*TODO*/][2 * i] *
-                                    strides[v];
-                              unsigned int ind2[VectorizedArrayType::size()];
-                              DEAL_II_OPENMP_SIMD_PRAGMA
-                              for (unsigned int v = 0;
-                                   v < VectorizedArrayType::size();
-                                   ++v)
-                                ind2[v] =
-                                  indices[v] +
-                                  index_array_hermite[0 /*TODO*/][2 * i + 1] *
-                                    strides[v];
-                              proc.hermite_grad_vectorized_indexed(
-                                temp1[i_],
-                                temp1[i_ + dofs_per_face],
-                                global_vector_ptr,
-                                global_vector_ptr,
-                                grad_weight,
-                                ind1,
-                                ind2);
-                            }
-                          else
-                            {
-                              Assert(false, ExcNotImplemented());
-                            }
-                        }
-                    else
+                    if (n_face_orientations == 1)
                       {
-                        if (integrate == false)
-                          for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
-                            temp1[i] = VectorizedArrayType();
-
-                        for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                          for (unsigned int i = 0; i < dofs_per_face; ++i)
-                            {
-                              const unsigned int i_ =
-                                reorientate(n_face_orientations == 1 ? 0 : v,
-                                            i);
-                              proc.hermite_grad(
-                                temp1[i_][v],
-                                temp1[i_ + dofs_per_face][v],
-                                global_vector_ptr
-                                  [indices[v] +
-                                   index_array_hermite
-                                       [n_face_orientations == 1 ? 0 : v]
-                                       [2 * i] *
-                                     strides[v]],
-                                global_vector_ptr
-                                  [indices[v] +
-                                   index_array_hermite
-                                       [n_face_orientations == 1 ? 0 : v]
-                                       [2 * i + 1] *
-                                     strides[v]],
-                                grad_weight[n_face_orientations == 1 ? 0 : v]);
-                            }
+                        const unsigned int i_ = reorientate(0, i);
+                        const unsigned int ind1 =
+                          index_array_hermite[0][2 * i] *
+                          VectorizedArrayType::size();
+                        const unsigned int ind2 =
+                          index_array_hermite[0][2 * i + 1] *
+                          VectorizedArrayType::size();
+                        proc.hermite_grad_vectorized_indexed(
+                          temp1[i_],
+                          temp1[i_ + dofs_per_face],
+                          vector_ptr + ind1,
+                          vector_ptr + ind2,
+                          grad_weight,
+                          indices,
+                          indices);
                       }
-                  }
-                else
-                  {
-                    if (n_filled_lanes == VectorizedArrayType::size())
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
-                        {
-                          if (n_face_orientations == 1)
-                            {
-                              unsigned int ind[VectorizedArrayType::size()];
-                              DEAL_II_OPENMP_SIMD_PRAGMA
-                              for (unsigned int v = 0;
-                                   v < VectorizedArrayType::size();
-                                   ++v)
-                                ind[v] = indices[v] +
-                                         index_array_nodal[0][i] * strides[v];
-                              const unsigned int i_ = reorientate(0, i);
-                              proc.value_vectorized_indexed(temp1[i_],
-                                                            global_vector_ptr,
-                                                            ind);
-                            }
-                          else
-                            {
-                              Assert(false, ExcNotImplemented());
-                            }
-                        }
                     else
                       {
-                        if (integrate == false)
-                          for (unsigned int i = 0; i < dofs_per_face; ++i)
-                            temp1[i] = VectorizedArrayType();
-
-                        for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                          for (unsigned int i = 0; i < dofs_per_face; ++i)
-                            proc.value(
-                              temp1[reorientate(
-                                n_face_orientations == 1 ? 0 : v, i)][v],
-                              global_vector_ptr
-                                [indices[v] +
-                                 index_array_nodal
-                                     [n_face_orientations == 1 ? 0 : v][i] *
-                                   strides[v]]);
+                        Assert(false, ExcNotImplemented());
                       }
                   }
               }
-
-            // case 4: contiguous indices without interleaving
-            else if ((n_face_orientations > 1 ||
-                      dof_info.index_storage_variants[dof_access_index][cell] ==
-                        MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                          contiguous))
+            else
               {
-                const unsigned int *indices =
-                  &dof_info.dof_indices_contiguous[dof_access_index]
-                                                  [cell *
-                                                   VectorizedArrayType::size()];
-                Number2_ *vector_ptr =
-                  global_vector_ptr + comp * static_dofs_per_component +
-                  dof_info
-                    .component_dof_indices_offset[active_fe_index]
-                                                 [first_selected_component];
-
-                const unsigned int n_filled_lanes =
-                  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
-
-                const bool vectorization_possible =
-                  (n_face_orientations == 1) &&
-                  (n_filled_lanes == VectorizedArrayType::size()) &&
-                  (sm_ptr != nullptr);
-
-                std::array<Number2_ *, VectorizedArrayType::size()>
-                  vector_ptrs = {};
-
-                if (vectorization_possible == false)
+                for (unsigned int i = 0; i < dofs_per_face; ++i)
                   {
                     if (n_face_orientations == 1)
                       {
-                        for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                          if (sm_ptr == nullptr)
-                            {
-                              vector_ptrs[v] = vector_ptr + indices[v];
-                            }
-                          else
-                            {
-                              const auto &temp =
-                                dof_info.dof_indices_contiguous_sm
-                                  [dof_access_index]
-                                  [cell * VectorizedArrayType::size() + v];
-                              vector_ptrs[v] = const_cast<Number *>(
-                                sm_ptr->operator[](temp.first).data() +
-                                temp.second + comp * static_dofs_per_component +
-                                dof_info.component_dof_indices_offset
-                                  [active_fe_index][first_selected_component]);
-                            }
-                      }
-                    else if (n_face_orientations == VectorizedArrayType::size())
-                      {
-                        for (unsigned int v = 0;
-                             v < VectorizedArrayType::size();
-                             ++v)
-                          if (cells[v] != numbers::invalid_unsigned_int)
-                            {
-                              if (sm_ptr == nullptr)
-                                {
-                                  vector_ptrs[v] =
-                                    vector_ptr +
-                                    dof_info
-                                      .dof_indices_contiguous[dof_access_index]
-                                                             [cells[v]];
-                                }
-                              else
-                                {
-                                  const auto &temp =
-                                    dof_info.dof_indices_contiguous_sm
-                                      [dof_access_index][cells[v]];
-                                  vector_ptrs[v] = const_cast<Number *>(
-                                    sm_ptr->operator[](temp.first).data() +
-                                    temp.second +
-                                    comp * static_dofs_per_component +
-                                    dof_info.component_dof_indices_offset
-                                      [active_fe_index]
-                                      [first_selected_component]);
-                                }
-                            }
+                        const unsigned int i_ = reorientate(0, i);
+                        const unsigned int ind =
+                          index_array_nodal[0][i] * VectorizedArrayType::size();
+                        proc.value_vectorized_indexed(temp1[i_],
+                                                      vector_ptr + ind,
+                                                      indices);
                       }
                     else
                       {
                         Assert(false, ExcNotImplemented());
                       }
                   }
+              }
+          }
 
-                if (do_gradients == true &&
-                    data.element_type ==
-                      MatrixFreeFunctions::tensor_symmetric_hermite)
-                  {
-                    if (vectorization_possible)
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+        // case 3: contiguous and interleaved indices with mixed stride
+        else if (n_face_orientations == 1 &&
+                 dof_info.index_storage_variants[dof_access_index][cell] ==
+                   MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                     interleaved_contiguous_mixed_strides)
+          {
+            AssertDimension(n_face_orientations, 1);
+
+            const unsigned int *strides =
+              &dof_info.dof_indices_interleave_strides
+                 [dof_access_index][cell * VectorizedArrayType::size()];
+            unsigned int indices[VectorizedArrayType::size()];
+            for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+              indices[v] =
+                dof_info.dof_indices_contiguous
+                  [dof_access_index][cell * VectorizedArrayType::size() + v] +
+                (dof_info
+                   .component_dof_indices_offset[active_fe_index]
+                                                [first_selected_component] +
+                 comp * static_dofs_per_component) *
+                  strides[v];
+            const unsigned int n_filled_lanes =
+              dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
+
+            if (fe_degree > 1 && do_gradients == true)
+              {
+                if (n_filled_lanes == VectorizedArrayType::size())
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      if (n_face_orientations == 1)
                         {
-                          const unsigned int ind1 =
-                            index_array_hermite[0][2 * i];
-                          const unsigned int ind2 =
-                            index_array_hermite[0][2 * i + 1];
                           const unsigned int i_ = reorientate(0, i);
-
+                          unsigned int       ind1[VectorizedArrayType::size()];
+                          DEAL_II_OPENMP_SIMD_PRAGMA
+                          for (unsigned int v = 0;
+                               v < VectorizedArrayType::size();
+                               ++v)
+                            ind1[v] =
+                              indices[v] +
+                              // TODO
+                              index_array_hermite[0][2 * i] * strides[v];
+                          unsigned int ind2[VectorizedArrayType::size()];
+                          DEAL_II_OPENMP_SIMD_PRAGMA
+                          for (unsigned int v = 0;
+                               v < VectorizedArrayType::size();
+                               ++v)
+                            ind2[v] =
+                              indices[v] +
+                              // TODO
+                              index_array_hermite[0][2 * i + 1] * strides[v];
                           proc.hermite_grad_vectorized_indexed(
                             temp1[i_],
                             temp1[i_ + dofs_per_face],
-                            vector_ptr + ind1,
-                            vector_ptr + ind2,
+                            global_vector_ptr,
+                            global_vector_ptr,
                             grad_weight,
-                            indices,
-                            indices);
+                            ind1,
+                            ind2);
                         }
-                    else if (n_face_orientations == 1)
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                      else
                         {
-                          const unsigned int ind1 =
-                            index_array_hermite[0][2 * i];
-                          const unsigned int ind2 =
-                            index_array_hermite[0][2 * i + 1];
-                          const unsigned int i_ = reorientate(0, i);
-
-                          for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                            proc.hermite_grad(temp1[i_][v],
-                                              temp1[i_ + dofs_per_face][v],
-                                              vector_ptrs[v][ind1],
-                                              vector_ptrs[v][ind2],
-                                              grad_weight[v]);
-
-                          if (integrate == false)
-                            for (unsigned int v = n_filled_lanes;
-                                 v < VectorizedArrayType::size();
-                                 ++v)
-                              {
-                                temp1[i_][v]                 = 0.0;
-                                temp1[i_ + dofs_per_face][v] = 0.0;
-                              }
+                          Assert(false, ExcNotImplemented());
                         }
-                    else
-                      {
-                        for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                          for (unsigned int i = 0; i < dofs_per_face; ++i)
-                            proc.hermite_grad(
-                              temp1[reorientate(v, i)][v],
-                              temp1[reorientate(v, i) + dofs_per_face][v],
-                              vector_ptrs[v][index_array_hermite[v][2 * i]],
-                              vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
-                              grad_weight[v]);
-                      }
-                  }
+                    }
                 else
                   {
-                    if (vectorization_possible)
+                    if (integrate == false)
+                      for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
+                        temp1[i] = VectorizedArrayType();
+
+                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
                       for (unsigned int i = 0; i < dofs_per_face; ++i)
                         {
-                          const unsigned int ind = index_array_nodal[0][i];
-                          const unsigned int i_  = reorientate(0, i);
-
+                          const unsigned int i_ =
+                            reorientate(n_face_orientations == 1 ? 0 : v, i);
+                          proc.hermite_grad(
+                            temp1[i_][v],
+                            temp1[i_ + dofs_per_face][v],
+                            global_vector_ptr
+                              [indices[v] +
+                               index_array_hermite
+                                   [n_face_orientations == 1 ? 0 : v][2 * i] *
+                                 strides[v]],
+                            global_vector_ptr
+                              [indices[v] +
+                               index_array_hermite[n_face_orientations == 1 ?
+                                                     0 :
+                                                     v][2 * i + 1] *
+                                 strides[v]],
+                            grad_weight[n_face_orientations == 1 ? 0 : v]);
+                        }
+                  }
+              }
+            else
+              {
+                if (n_filled_lanes == VectorizedArrayType::size())
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      if (n_face_orientations == 1)
+                        {
+                          unsigned int ind[VectorizedArrayType::size()];
+                          DEAL_II_OPENMP_SIMD_PRAGMA
+                          for (unsigned int v = 0;
+                               v < VectorizedArrayType::size();
+                               ++v)
+                            ind[v] =
+                              indices[v] + index_array_nodal[0][i] * strides[v];
+                          const unsigned int i_ = reorientate(0, i);
                           proc.value_vectorized_indexed(temp1[i_],
-                                                        vector_ptr + ind,
-                                                        indices);
+                                                        global_vector_ptr,
+                                                        ind);
                         }
-                    else if (n_face_orientations == 1)
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                      else
                         {
-                          const unsigned int ind = index_array_nodal[0][i];
-                          const unsigned int i_  = reorientate(0, i);
+                          Assert(false, ExcNotImplemented());
+                        }
+                    }
+                else
+                  {
+                    if (integrate == false)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        temp1[i] = VectorizedArrayType();
 
-                          for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                            proc.value(temp1[i_][v], vector_ptrs[v][ind]);
+                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        proc.value(
+                          temp1[reorientate(n_face_orientations == 1 ? 0 : v,
+                                            i)][v],
+                          global_vector_ptr
+                            [indices[v] +
+                             index_array_nodal[n_face_orientations == 1 ? 0 : v]
+                                              [i] *
+                               strides[v]]);
+                  }
+              }
+          }
 
-                          if (integrate == false)
-                            for (unsigned int v = n_filled_lanes;
-                                 v < VectorizedArrayType::size();
-                                 ++v)
-                              temp1[i_][v] = 0.0;
+        // case 4: contiguous indices without interleaving
+        else if ((n_face_orientations > 1 ||
+                  dof_info.index_storage_variants[dof_access_index][cell] ==
+                    MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                      contiguous))
+          {
+            const unsigned int *indices =
+              &dof_info
+                 .dof_indices_contiguous[dof_access_index]
+                                        [cell * VectorizedArrayType::size()];
+            Number2_ *vector_ptr =
+              global_vector_ptr + comp * static_dofs_per_component +
+              dof_info.component_dof_indices_offset[active_fe_index]
+                                                   [first_selected_component];
+
+            const unsigned int n_filled_lanes =
+              dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
+
+            const bool vectorization_possible =
+              (n_face_orientations == 1) &&
+              (n_filled_lanes == VectorizedArrayType::size()) &&
+              (sm_ptr != nullptr);
+
+            std::array<Number2_ *, VectorizedArrayType::size()> vector_ptrs =
+              {};
+
+            if (vectorization_possible == false)
+              {
+                if (n_face_orientations == 1)
+                  {
+                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                      if (sm_ptr == nullptr)
+                        {
+                          vector_ptrs[v] = vector_ptr + indices[v];
                         }
-                    else
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                      else
                         {
-                          for (unsigned int v = 0;
-                               v < VectorizedArrayType::size();
-                               ++v)
-                            if (cells[v] != numbers::invalid_unsigned_int)
-                              proc.value(
-                                temp1[reorientate(v, i)][v],
-                                vector_ptrs[v][index_array_nodal[v][i]]);
-                            else if (integrate == false)
-                              temp1[i][v] = 0.0;
+                          const auto &temp =
+                            dof_info.dof_indices_contiguous_sm
+                              [dof_access_index]
+                              [cell * VectorizedArrayType::size() + v];
+                          vector_ptrs[v] = const_cast<Number *>(
+                            sm_ptr->operator[](temp.first).data() +
+                            temp.second + comp * static_dofs_per_component +
+                            dof_info.component_dof_indices_offset
+                              [active_fe_index][first_selected_component]);
+                        }
+                  }
+                else if (n_face_orientations == VectorizedArrayType::size())
+                  {
+                    for (unsigned int v = 0; v < VectorizedArrayType::size();
+                         ++v)
+                      if (cells[v] != numbers::invalid_unsigned_int)
+                        {
+                          if (sm_ptr == nullptr)
+                            {
+                              vector_ptrs[v] =
+                                vector_ptr +
+                                dof_info
+                                  .dof_indices_contiguous[dof_access_index]
+                                                         [cells[v]];
+                            }
+                          else
+                            {
+                              const auto &temp =
+                                dof_info
+                                  .dof_indices_contiguous_sm[dof_access_index]
+                                                            [cells[v]];
+                              vector_ptrs[v] = const_cast<Number *>(
+                                sm_ptr->operator[](temp.first).data() +
+                                temp.second + comp * static_dofs_per_component +
+                                dof_info.component_dof_indices_offset
+                                  [active_fe_index][first_selected_component]);
+                            }
                         }
                   }
+                else
+                  {
+                    Assert(false, ExcNotImplemented());
+                  }
+              }
+
+            if (do_gradients == true &&
+                data.element_type ==
+                  MatrixFreeFunctions::tensor_symmetric_hermite)
+              {
+                if (vectorization_possible)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      const unsigned int ind1 = index_array_hermite[0][2 * i];
+                      const unsigned int ind2 =
+                        index_array_hermite[0][2 * i + 1];
+                      const unsigned int i_ = reorientate(0, i);
+
+                      proc.hermite_grad_vectorized_indexed(
+                        temp1[i_],
+                        temp1[i_ + dofs_per_face],
+                        vector_ptr + ind1,
+                        vector_ptr + ind2,
+                        grad_weight,
+                        indices,
+                        indices);
+                    }
+                else if (n_face_orientations == 1)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      const unsigned int ind1 = index_array_hermite[0][2 * i];
+                      const unsigned int ind2 =
+                        index_array_hermite[0][2 * i + 1];
+                      const unsigned int i_ = reorientate(0, i);
+
+                      for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                        proc.hermite_grad(temp1[i_][v],
+                                          temp1[i_ + dofs_per_face][v],
+                                          vector_ptrs[v][ind1],
+                                          vector_ptrs[v][ind2],
+                                          grad_weight[v]);
+
+                      if (integrate == false)
+                        for (unsigned int v = n_filled_lanes;
+                             v < VectorizedArrayType::size();
+                             ++v)
+                          {
+                            temp1[i_][v]                 = 0.0;
+                            temp1[i_ + dofs_per_face][v] = 0.0;
+                          }
+                    }
+                else
+                  {
+                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        proc.hermite_grad(
+                          temp1[reorientate(v, i)][v],
+                          temp1[reorientate(v, i) + dofs_per_face][v],
+                          vector_ptrs[v][index_array_hermite[v][2 * i]],
+                          vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
+                          grad_weight[v]);
+                  }
               }
             else
               {
-                // case 5: default vector access
-                // for the integrate_scatter path (integrate == true), we
-                // need to only prepare the data in this function for all
-                // components to later call distribute_local_to_global();
-                // for the gather_evaluate path (integrate == false), we
-                // instead want to leave early because we need to get the
-                // vector data from somewhere else
-                proc.default_operation(temp1, comp);
-                if (integrate)
-                  accesses_global_vector = false;
+                if (vectorization_possible)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      const unsigned int ind = index_array_nodal[0][i];
+                      const unsigned int i_  = reorientate(0, i);
+
+                      proc.value_vectorized_indexed(temp1[i_],
+                                                    vector_ptr + ind,
+                                                    indices);
+                    }
+                else if (n_face_orientations == 1)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      const unsigned int ind = index_array_nodal[0][i];
+                      const unsigned int i_  = reorientate(0, i);
+
+                      for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                        proc.value(temp1[i_][v], vector_ptrs[v][ind]);
+
+                      if (integrate == false)
+                        for (unsigned int v = n_filled_lanes;
+                             v < VectorizedArrayType::size();
+                             ++v)
+                          temp1[i_][v] = 0.0;
+                    }
                 else
-                  return false;
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      for (unsigned int v = 0; v < VectorizedArrayType::size();
+                           ++v)
+                        if (cells[v] != numbers::invalid_unsigned_int)
+                          proc.value(temp1[reorientate(v, i)][v],
+                                     vector_ptrs[v][index_array_nodal[v][i]]);
+                    }
               }
           }
         else
           {
-            // case 5: default vector access
-            proc.default_operation(temp1, comp);
-            if (integrate)
-              accesses_global_vector = false;
-            else
-              return false;
+            // We should not end up here, this should be caught by
+            // FEFaceEvaluationImplGatherEvaluateSelector::supports()
+            Assert(false, ExcInternalError());
           }
-
-        if (!integrate)
-          proc.in_face_operation(temp1, comp);
       }
-
-    if (!integrate &&
-        (face_orientations[0] > 0 &&
-         subface_index < GeometryInfo<dim>::max_children_per_cell))
-      {
-        AssertDimension(n_face_orientations, 1);
-        adjust_for_face_orientation(dim,
-                                    n_components,
-                                    face_orientations[0],
-                                    orientation_map,
-                                    false,
-                                    do_values,
-                                    do_gradients,
-                                    do_hessians,
-                                    data.n_q_points_face,
-                                    scratch_data,
-                                    values_quad,
-                                    gradients_quad,
-                                    hessians_quad);
-      }
-
-    return accesses_global_vector;
+    */
   }
 
 
@@ -4065,121 +3701,118 @@ namespace internal
     template <int fe_degree, int n_q_points_1d>
     static bool
     run(const unsigned int                          n_components,
-        const unsigned int                          n_face_orientations,
+        const EvaluationFlags::EvaluationFlags      evaluation_flag,
         const Number2 *                             src_ptr,
         const std::vector<ArrayView<const Number>> *sm_ptr,
-        const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-        const MatrixFreeFunctions::DoFInfo &                       dof_info,
-        VectorizedArrayType *                                      values_quad,
-        VectorizedArrayType *gradients_quad,
-        VectorizedArrayType *hessians_quad,
-        VectorizedArrayType *scratch_data,
-        const bool           evaluate_values,
-        const bool           evaluate_gradients,
-        const bool           evaluate_hessians,
-        const unsigned int   active_fe_index,
-        const unsigned int   first_selected_component,
-        const std::array<unsigned int, VectorizedArrayType::size()> cells,
-        const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
-        const unsigned int                                 subface_index,
-        const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
-        const std::array<unsigned int, VectorizedArrayType::size()>
-                                      face_orientations,
-        const Table<2, unsigned int> &orientation_map)
+        FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval)
     {
-      if (src_ptr == nullptr)
-        return false;
+      Assert(fe_degree > -1, ExcInternalError());
+      Assert(eval.get_shape_info().element_type <=
+               MatrixFreeFunctions::tensor_symmetric,
+             ExcInternalError());
 
-      if (data.element_type == MatrixFreeFunctions::tensor_none)
-        return false;
+      const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
+
+      VectorizedArrayType *temp = eval.get_scratch_data().begin();
+      VectorizedArrayType *scratch_data =
+        temp + 3 * n_components * dofs_per_face;
+
+      Processor<fe_degree, n_q_points_1d> p;
+
+      if (eval.get_dof_access_index() ==
+            internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+          eval.get_is_interior_face() == false)
+        fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
+          p, n_components, evaluation_flag, src_ptr, sm_ptr, eval, temp);
+      else
+        fe_face_evaluation_process_and_io<1>(
+          p, n_components, evaluation_flag, src_ptr, sm_ptr, eval, temp);
+
+      const unsigned int subface_index = eval.get_subface_index();
+
+      if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
+        FEFaceEvaluationImpl<
+          true,
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          VectorizedArrayType>::evaluate_in_face(n_components,
+                                                 evaluation_flag,
+                                                 eval.get_shape_info().data[0],
+                                                 temp,
+                                                 eval.begin_values(),
+                                                 eval.begin_gradients(),
+                                                 eval.begin_hessians(),
+                                                 scratch_data,
+                                                 subface_index);
+      else
+        FEFaceEvaluationImpl<
+          false,
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          VectorizedArrayType>::evaluate_in_face(n_components,
+                                                 evaluation_flag,
+                                                 eval.get_shape_info().data[0],
+                                                 temp,
+                                                 eval.begin_values(),
+                                                 eval.begin_gradients(),
+                                                 eval.begin_hessians(),
+                                                 scratch_data,
+                                                 subface_index);
+
+      if (eval.get_face_orientation())
+        adjust_for_face_orientation(
+          dim,
+          n_components,
+          evaluation_flag,
+          &eval.get_orientation_map()(eval.get_face_orientation(), 0),
+          false,
+          eval.get_shape_info().n_q_points_face,
+          temp,
+          eval.begin_values(),
+          eval.begin_gradients(),
+          eval.begin_hessians());
+
+      return false;
+    }
 
-      (void)sm_ptr;
-
-      Processor<fe_degree, n_q_points_1d> p(n_components,
-                                            false,
-                                            src_ptr,
-                                            sm_ptr,
-                                            data,
-                                            dof_info,
-                                            values_quad,
-                                            gradients_quad,
-                                            hessians_quad,
-                                            scratch_data,
-                                            evaluate_values,
-                                            evaluate_gradients,
-                                            evaluate_hessians,
-                                            active_fe_index,
-                                            first_selected_component,
-                                            cells,
-                                            face_nos,
-                                            subface_index,
-                                            dof_access_index,
-                                            face_orientations,
-                                            orientation_map);
-
-      if (n_face_orientations == VectorizedArrayType::size())
-        return fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
-          p);
+    static bool
+    supports(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
+        &                                                          shape_info,
+      const Number *                                               vector_ptr,
+      internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
+    {
+      const unsigned int fe_degree = shape_info.data[0].fe_degree;
+      if (fe_degree < 1 || !shape_info.data[0].nodal_at_cell_boundaries ||
+          (evaluation_flag & EvaluationFlags::gradients &&
+           (fe_degree < 2 ||
+            shape_info.data[0].element_type !=
+              MatrixFreeFunctions::tensor_symmetric_hermite)) ||
+          (evaluation_flag & EvaluationFlags::hessians) ||
+          vector_ptr == nullptr ||
+          shape_info.data[0].element_type >
+            MatrixFreeFunctions::tensor_symmetric ||
+          storage <
+            MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
+        return false;
       else
-        return fe_face_evaluation_process_and_io<1>(p);
+        return false;
     }
 
   private:
     template <int fe_degree, int n_q_points_1d>
     struct Processor
     {
-      static const int dim_           = dim;
-      static const int fe_degree_     = fe_degree;
-      static const int n_q_points_1d_ = n_q_points_1d;
-      using VectorizedArrayType_      = VectorizedArrayType;
-      using Number_                   = Number;
-      using Number2_                  = const Number2;
-
-      Processor(
-        const unsigned int                          n_components,
-        const bool                                  integrate,
-        const Number2 *                             global_vector_ptr,
-        const std::vector<ArrayView<const Number>> *sm_ptr,
-        const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-        const MatrixFreeFunctions::DoFInfo &                       dof_info,
-        VectorizedArrayType *                                      values_quad,
-        VectorizedArrayType *gradients_quad,
-        VectorizedArrayType *hessians_quad,
-        VectorizedArrayType *scratch_data,
-        const bool           do_values,
-        const bool           do_gradients,
-        const bool           do_hessians,
-        const unsigned int   active_fe_index,
-        const unsigned int   first_selected_component,
-        const std::array<unsigned int, VectorizedArrayType::size()> cells,
-        const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
-        const unsigned int                                 subface_index,
-        const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
-        const std::array<unsigned int, VectorizedArrayType::size()>
-                                      face_orientations,
-        const Table<2, unsigned int> &orientation_map)
-        : n_components(n_components)
-        , integrate(integrate)
-        , global_vector_ptr(global_vector_ptr)
-        , sm_ptr(sm_ptr)
-        , data(data)
-        , dof_info(dof_info)
-        , values_quad(values_quad)
-        , gradients_quad(gradients_quad)
-        , hessians_quad(hessians_quad)
-        , scratch_data(scratch_data)
-        , do_values(do_values)
-        , do_gradients(do_gradients)
-        , do_hessians(do_hessians)
-        , active_fe_index(active_fe_index)
-        , first_selected_component(first_selected_component)
-        , cells(cells)
-        , face_nos(face_nos)
-        , subface_index(subface_index)
-        , dof_access_index(dof_access_index)
-        , face_orientations(face_orientations)
-        , orientation_map(orientation_map)
-      {}
+      static const bool do_integrate   = false;
+      static const int  dim_           = dim;
+      static const int  fe_degree_     = fe_degree;
+      static const int  n_q_points_1d_ = n_q_points_1d;
+      using VectorizedArrayType_       = VectorizedArrayType;
+      using Number_                    = Number;
+      using Number2_                   = const Number2;
 
       template <typename T0, typename T1, typename T2>
       void
@@ -4243,86 +3876,6 @@ namespace internal
         // case 3b)
         temp = src_ptr;
       }
-
-      template <typename T1>
-      void
-      default_operation(const T1 &, const unsigned int)
-      {
-        // case 5)
-      }
-
-      template <typename T1>
-      void
-      in_face_operation(T1 &temp1, const unsigned int comp)
-      {
-        const unsigned int dofs_per_face =
-          fe_degree > -1 ?
-            Utilities::pow(fe_degree + 1, dim - 1) :
-            Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
-        const unsigned int n_q_points =
-          fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
-                           data.n_q_points_face;
-        constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
-        if (fe_degree > -1 &&
-            subface_index >= GeometryInfo<dim>::max_children_per_cell &&
-            data.element_type <= MatrixFreeFunctions::tensor_symmetric)
-          FEFaceEvaluationImpl<true,
-                               dim,
-                               fe_degree,
-                               n_q_points_1d,
-                               VectorizedArrayType>::
-            evaluate_in_face(/* n_components */ 1,
-                             data,
-                             temp1,
-                             values_quad + comp * n_q_points,
-                             gradients_quad + comp * dim * n_q_points,
-                             hessians_quad + comp * hdim * n_q_points,
-                             scratch_data + 3 * dofs_per_face,
-                             do_values,
-                             do_gradients,
-                             do_hessians,
-                             subface_index);
-        else
-          FEFaceEvaluationImpl<false,
-                               dim,
-                               fe_degree,
-                               n_q_points_1d,
-                               VectorizedArrayType>::
-            evaluate_in_face(/* n_components */ 1,
-                             data,
-                             temp1,
-                             values_quad + comp * n_q_points,
-                             gradients_quad + comp * dim * n_q_points,
-                             hessians_quad + comp * hdim * n_q_points,
-                             scratch_data + 3 * dofs_per_face,
-                             do_values,
-                             do_gradients,
-                             do_hessians,
-                             subface_index);
-      }
-
-      const unsigned int                          n_components;
-      const bool                                  integrate;
-      const Number2 *                             global_vector_ptr;
-      const std::vector<ArrayView<const Number>> *sm_ptr;
-      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data;
-      const MatrixFreeFunctions::DoFInfo &                       dof_info;
-      VectorizedArrayType *                                      values_quad;
-      VectorizedArrayType *                                      gradients_quad;
-      VectorizedArrayType *                                      hessians_quad;
-      VectorizedArrayType *                                      scratch_data;
-      const bool                                                 do_values;
-      const bool                                                 do_gradients;
-      const bool                                                 do_hessians;
-      const unsigned int active_fe_index;
-      const unsigned int first_selected_component;
-      const std::array<unsigned int, VectorizedArrayType::size()> cells;
-      const std::array<unsigned int, VectorizedArrayType::size()> face_nos;
-      const unsigned int                                          subface_index;
-      const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index;
-      const std::array<unsigned int, VectorizedArrayType::size()>
-                                    face_orientations;
-      const Table<2, unsigned int> &orientation_map;
     };
   };
 
@@ -4335,147 +3888,94 @@ namespace internal
     template <int fe_degree, int n_q_points_1d>
     static bool
     run(const unsigned int                           n_components,
-        const unsigned int                           n_face_orientations,
+        const EvaluationFlags::EvaluationFlags       integration_flag,
         Number2 *                                    dst_ptr,
         const std::vector<ArrayView<const Number2>> *sm_ptr,
-        const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-        const MatrixFreeFunctions::DoFInfo &                       dof_info,
-        VectorizedArrayType *                                      values_array,
-        VectorizedArrayType *                                      values_quad,
-        VectorizedArrayType *gradients_quad,
-        VectorizedArrayType *hessians_quad,
-        VectorizedArrayType *scratch_data,
-        const bool           integrate_values,
-        const bool           integrate_gradients,
-        const bool           integrate_hessians,
-        const unsigned int   active_fe_index,
-        const unsigned int   first_selected_component,
-        const std::array<unsigned int, VectorizedArrayType::size()> cells,
-        const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
-        const unsigned int                                 subface_index,
-        const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
-        const std::array<unsigned int, VectorizedArrayType::size()>
-                                      face_orientations,
-        const Table<2, unsigned int> &orientation_map)
+        FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval)
     {
-      (void)sm_ptr;
+      Assert(fe_degree > -1, ExcInternalError());
+      Assert(eval.get_shape_info().element_type <=
+               MatrixFreeFunctions::tensor_symmetric,
+             ExcInternalError());
 
-      if (dst_ptr == nullptr ||
-          data.element_type == MatrixFreeFunctions::tensor_none)
-        {
-          AssertDimension(n_face_orientations, 1);
-
-          // for block vectors simply integrate
-          FEFaceEvaluationImplIntegrateSelector<dim, VectorizedArrayType>::
-            template run<fe_degree, n_q_points_1d>(n_components,
-                                                   data,
-                                                   values_array,
-                                                   values_quad,
-                                                   gradients_quad,
-                                                   hessians_quad,
-                                                   scratch_data,
-                                                   integrate_values,
-                                                   integrate_gradients,
-                                                   integrate_hessians,
-                                                   face_nos[0],
-                                                   subface_index,
-                                                   face_orientations[0],
-                                                   orientation_map);
-
-          // default vector access
-          return false;
-        }
+      const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
+
+      VectorizedArrayType *temp = eval.get_scratch_data().begin();
+      VectorizedArrayType *scratch_data =
+        temp + 3 * n_components * dofs_per_face;
+
+      if (eval.get_face_orientation())
+        adjust_for_face_orientation(
+          dim,
+          n_components,
+          integration_flag,
+          &eval.get_orientation_map()(eval.get_face_orientation(), 0),
+          true,
+          eval.get_shape_info().n_q_points_face,
+          temp,
+          eval.begin_values(),
+          eval.begin_gradients(),
+          eval.begin_hessians());
+
+      const unsigned int subface_index = eval.get_subface_index();
+
+      if (fe_degree > -1 && eval.get_subface_index() >=
+                              GeometryInfo<dim - 1>::max_children_per_cell)
+        FEFaceEvaluationImpl<
+          true,
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          VectorizedArrayType>::integrate_in_face(n_components,
+                                                  integration_flag,
+                                                  eval.get_shape_info().data[0],
+                                                  temp,
+                                                  eval.begin_values(),
+                                                  eval.begin_gradients(),
+                                                  eval.begin_hessians(),
+                                                  scratch_data,
+                                                  subface_index);
+      else
+        FEFaceEvaluationImpl<
+          false,
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          VectorizedArrayType>::integrate_in_face(n_components,
+                                                  integration_flag,
+                                                  eval.get_shape_info().data[0],
+                                                  temp,
+                                                  eval.begin_values(),
+                                                  eval.begin_gradients(),
+                                                  eval.begin_hessians(),
+                                                  scratch_data,
+                                                  subface_index);
 
+      Processor<fe_degree, n_q_points_1d> p;
 
-      Processor<fe_degree, n_q_points_1d> p(values_array,
-                                            n_components,
-                                            true,
-                                            dst_ptr,
-                                            sm_ptr,
-                                            data,
-                                            dof_info,
-                                            values_quad,
-                                            gradients_quad,
-                                            hessians_quad,
-                                            scratch_data,
-                                            integrate_values,
-                                            integrate_gradients,
-                                            integrate_hessians,
-                                            active_fe_index,
-                                            first_selected_component,
-                                            cells,
-                                            face_nos,
-                                            subface_index,
-                                            dof_access_index,
-                                            face_orientations,
-                                            orientation_map);
-
-      if (n_face_orientations == VectorizedArrayType::size())
-        return fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
-          p);
+      if (eval.get_dof_access_index() ==
+            internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+          eval.get_is_interior_face() == false)
+        fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
+          p, n_components, integration_flag, dst_ptr, sm_ptr, eval, temp);
       else
-        return fe_face_evaluation_process_and_io<1>(p);
+        fe_face_evaluation_process_and_io<1>(
+          p, n_components, integration_flag, dst_ptr, sm_ptr, eval, temp);
+
+      return false;
     }
 
   private:
     template <int fe_degree, int n_q_points_1d>
     struct Processor
     {
-      static const int dim_           = dim;
-      static const int fe_degree_     = fe_degree;
-      static const int n_q_points_1d_ = n_q_points_1d;
-      using VectorizedArrayType_      = VectorizedArrayType;
-      using Number_                   = Number;
-      using Number2_                  = Number2;
-
-
-      Processor(
-        VectorizedArrayType *                       values_array,
-        const unsigned int                          n_components,
-        const bool                                  integrate,
-        Number2 *                                   global_vector_ptr,
-        const std::vector<ArrayView<const Number>> *sm_ptr,
-        const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-        const MatrixFreeFunctions::DoFInfo &                       dof_info,
-        VectorizedArrayType *                                      values_quad,
-        VectorizedArrayType *gradients_quad,
-        VectorizedArrayType *hessians_quad,
-        VectorizedArrayType *scratch_data,
-        const bool           do_values,
-        const bool           do_gradients,
-        const bool           do_hessians,
-        const unsigned int   active_fe_index,
-        const unsigned int   first_selected_component,
-        const std::array<unsigned int, VectorizedArrayType::size()> cells,
-        const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
-        const unsigned int                                 subface_index,
-        const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
-        const std::array<unsigned int, VectorizedArrayType::size()>
-                                      face_orientations,
-        const Table<2, unsigned int> &orientation_map)
-        : values_array(values_array)
-        , n_components(n_components)
-        , integrate(integrate)
-        , global_vector_ptr(global_vector_ptr)
-        , sm_ptr(sm_ptr)
-        , data(data)
-        , dof_info(dof_info)
-        , values_quad(values_quad)
-        , gradients_quad(gradients_quad)
-        , hessians_quad(hessians_quad)
-        , scratch_data(scratch_data)
-        , do_values(do_values)
-        , do_gradients(do_gradients)
-        , do_hessians(do_hessians)
-        , active_fe_index(active_fe_index)
-        , first_selected_component(first_selected_component)
-        , cells(cells)
-        , face_nos(face_nos)
-        , subface_index(subface_index)
-        , dof_access_index(dof_access_index)
-        , face_orientations(face_orientations)
-        , orientation_map(orientation_map)
-      {}
+      static const bool do_integrate   = true;
+      static const int  dim_           = dim;
+      static const int  fe_degree_     = fe_degree;
+      static const int  n_q_points_1d_ = n_q_points_1d;
+      using VectorizedArrayType_       = VectorizedArrayType;
+      using Number_                    = Number;
+      using Number2_                   = Number2;
 
       template <typename T0, typename T1, typename T2, typename T3, typename T4>
       void
@@ -4547,101 +4047,6 @@ namespace internal
         // case 3b)
         dst_ptr += temp;
       }
-
-      template <typename T0>
-      void
-      default_operation(const T0 &temp1, const unsigned int comp)
-      {
-        // case 5: default vector access, must be handled separately, just do
-        // the face-normal interpolation
-
-        FEFaceNormalEvaluationImpl<dim, fe_degree, VectorizedArrayType>::
-          template interpolate<false, false>(
-            /* n_components */ 1,
-            data,
-            temp1,
-            values_array + comp * data.dofs_per_component_on_cell,
-            do_gradients,
-            do_hessians,
-            face_nos[0]);
-      }
-
-      template <typename T0>
-      void
-      in_face_operation(T0 &temp1, const unsigned int comp)
-      {
-        const unsigned int dofs_per_face =
-          fe_degree > -1 ?
-            Utilities::pow(fe_degree + 1, dim - 1) :
-            Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
-        const unsigned int n_q_points =
-          fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
-                           data.n_q_points_face;
-        constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
-        if (fe_degree > -1 &&
-            subface_index >= GeometryInfo<dim>::max_children_per_cell &&
-            data.element_type <=
-              internal::MatrixFreeFunctions::tensor_symmetric)
-          internal::FEFaceEvaluationImpl<true,
-                                         dim,
-                                         fe_degree,
-                                         n_q_points_1d,
-                                         VectorizedArrayType>::
-            integrate_in_face(/* n_components */ 1,
-                              data,
-                              temp1,
-                              values_quad + comp * n_q_points,
-                              gradients_quad + dim * comp * n_q_points,
-                              hessians_quad + hdim * comp * n_q_points,
-                              scratch_data + 3 * dofs_per_face,
-                              do_values,
-                              do_gradients,
-                              do_hessians,
-                              subface_index);
-        else
-          internal::FEFaceEvaluationImpl<false,
-                                         dim,
-                                         fe_degree,
-                                         n_q_points_1d,
-                                         VectorizedArrayType>::
-            integrate_in_face(/* n_components */ 1,
-                              data,
-                              temp1,
-                              values_quad + comp * n_q_points,
-                              gradients_quad + dim * comp * n_q_points,
-                              hessians_quad + hdim * comp * n_q_points,
-                              scratch_data + 3 * dofs_per_face,
-                              do_values,
-                              do_gradients,
-                              do_hessians,
-                              subface_index);
-      }
-
-      VectorizedArrayType *values_array;
-
-
-      const unsigned int                          n_components;
-      const bool                                  integrate;
-      Number2 *                                   global_vector_ptr;
-      const std::vector<ArrayView<const Number>> *sm_ptr;
-      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data;
-      const MatrixFreeFunctions::DoFInfo &                       dof_info;
-      VectorizedArrayType *                                      values_quad;
-      VectorizedArrayType *                                      gradients_quad;
-      VectorizedArrayType *                                      hessians_quad;
-      VectorizedArrayType *                                      scratch_data;
-      const bool                                                 do_values;
-      const bool                                                 do_gradients;
-      const bool                                                 do_hessians;
-      const unsigned int active_fe_index;
-      const unsigned int first_selected_component;
-      const std::array<unsigned int, VectorizedArrayType::size()> cells;
-      const std::array<unsigned int, VectorizedArrayType::size()> face_nos;
-      const unsigned int                                          subface_index;
-      const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index;
-      const std::array<unsigned int, VectorizedArrayType::size()>
-                                    face_orientations;
-      const Table<2, unsigned int> &orientation_map;
     };
   };
 
index 1c9d8b61efc9ae8e077ffefa7a9ebe46ea32cb2f..e46101578645c58f045762b4ae1c82d0a8faced0 100644 (file)
@@ -35,6 +35,9 @@ DEAL_II_NAMESPACE_OPEN
 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
 struct SelectEvaluator
 {
+  using ScalarNumber =
+    typename internal::VectorizedArrayTrait<Number>::value_type;
+
   /**
    * Chooses an appropriate evaluation strategy for the evaluate function, i.e.
    * this calls internal::FEEvaluationImpl::evaluate(),
@@ -45,12 +48,8 @@ struct SelectEvaluator
   static void
   evaluate(const unsigned int                     n_components,
            const EvaluationFlags::EvaluationFlags evaluation_flag,
-           const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-           Number *values_dofs_actual,
-           Number *values_quad,
-           Number *gradients_quad,
-           Number *hessians_quad,
-           Number *scratch_data);
+           const Number *                         values_dofs,
+           FEEvaluationBaseData<dim, ScalarNumber, false, Number> &eval);
 
   /**
    * Chooses an appropriate evaluation strategy for the integrate function, i.e.
@@ -62,12 +61,8 @@ struct SelectEvaluator
   static void
   integrate(const unsigned int                     n_components,
             const EvaluationFlags::EvaluationFlags integration_flag,
-            const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-            Number *   values_dofs_actual,
-            Number *   values_quad,
-            Number *   gradients_quad,
-            Number *   hessians_quad,
-            Number *   scratch_data,
+            Number *                               values_dofs,
+            FEEvaluationBaseData<dim, ScalarNumber, false, Number> &eval,
             const bool sum_into_values_array = false);
 };
 
@@ -79,24 +74,14 @@ inline void
 SelectEvaluator<dim, fe_degree, n_q_points_1d, Number>::evaluate(
   const unsigned int                                      n_components,
   const EvaluationFlags::EvaluationFlags                  evaluation_flag,
-  const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-  Number *                                                values_dofs_actual,
-  Number *                                                values_quad,
-  Number *                                                gradients_quad,
-  Number *                                                hessians_quad,
-  Number *                                                scratch_data)
+  const Number *                                          values_dofs,
+  FEEvaluationBaseData<dim, ScalarNumber, false, Number> &eval)
 {
   Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError());
 
-  internal::FEEvaluationImplEvaluateSelector<dim, Number>::
-    template run<fe_degree, n_q_points_1d>(n_components,
-                                           evaluation_flag,
-                                           shape_info,
-                                           values_dofs_actual,
-                                           values_quad,
-                                           gradients_quad,
-                                           hessians_quad,
-                                           scratch_data);
+  internal::FEEvaluationImplEvaluateSelector<dim, Number>::template run<
+    fe_degree,
+    n_q_points_1d>(n_components, evaluation_flag, values_dofs, eval);
 }
 
 
@@ -106,26 +91,15 @@ inline void
 SelectEvaluator<dim, fe_degree, n_q_points_1d, Number>::integrate(
   const unsigned int                                      n_components,
   const EvaluationFlags::EvaluationFlags                  integration_flag,
-  const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-  Number *                                                values_dofs_actual,
-  Number *                                                values_quad,
-  Number *                                                gradients_quad,
-  Number *                                                hessians_quad,
-  Number *                                                scratch_data,
+  Number *                                                values_dofs,
+  FEEvaluationBaseData<dim, ScalarNumber, false, Number> &eval,
   const bool                                              sum_into_values_array)
 {
   Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError());
 
   internal::FEEvaluationImplIntegrateSelector<dim, Number>::
-    template run<fe_degree, n_q_points_1d>(n_components,
-                                           integration_flag,
-                                           shape_info,
-                                           values_dofs_actual,
-                                           values_quad,
-                                           gradients_quad,
-                                           hessians_quad,
-                                           scratch_data,
-                                           sum_into_values_array);
+    template run<fe_degree, n_q_points_1d>(
+      n_components, integration_flag, values_dofs, eval, sum_into_values_array);
 }
 #endif // DOXYGEN
 
index 1dbeef01a360ae19c9ac48ee63dc305a50358171..e0152471f7268204beb512553bc9483b9557077a 100644 (file)
@@ -35,33 +35,23 @@ class FEEvaluationBaseData;
 
 namespace internal
 {
-  template <int dim,
-            typename Number,
-            typename VectorizedArrayType = VectorizedArray<Number>>
+  template <int dim, typename Number, typename VectorizedArrayType>
   struct FEEvaluationFactory
   {
     static void
     evaluate(
       const unsigned int                     n_components,
       const EvaluationFlags::EvaluationFlags evaluation_flag,
-      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
-      VectorizedArrayType *values_dofs_actual,
-      VectorizedArrayType *values_quad,
-      VectorizedArrayType *gradients_quad,
-      VectorizedArrayType *hessians_quad,
-      VectorizedArrayType *scratch_data);
+      const VectorizedArrayType *            values_dofs,
+      FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> &eval);
 
     static void
     integrate(
       const unsigned int                     n_components,
       const EvaluationFlags::EvaluationFlags integration_flag,
-      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
-      VectorizedArrayType *values_dofs_actual,
-      VectorizedArrayType *values_quad,
-      VectorizedArrayType *gradients_quad,
-      VectorizedArrayType *hessians_quad,
-      VectorizedArrayType *scratch_data,
-      const bool           sum_into_values_array);
+      VectorizedArrayType *                  values_dofs,
+      FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> &eval,
+      const bool sum_into_values_array);
 
     static bool
     fast_evaluation_supported(const unsigned int given_degree,
@@ -76,87 +66,34 @@ namespace internal
   struct FEFaceEvaluationFactory
   {
     static void
-    evaluate(const unsigned int n_components,
-             const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-             const VectorizedArrayType *   values_array,
-             VectorizedArrayType *         values_quad,
-             VectorizedArrayType *         gradients_quad,
-             VectorizedArrayType *         hessians_quad,
-             VectorizedArrayType *         scratch_data,
-             const bool                    evaluate_values,
-             const bool                    evaluate_gradients,
-             const bool                    evaluate_hessians,
-             const unsigned int            face_no,
-             const unsigned int            subface_index,
-             const unsigned int            face_orientation,
-             const Table<2, unsigned int> &orientation_map);
+    evaluate(
+      const unsigned int                     n_components,
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      const VectorizedArrayType *            values_dofs,
+      FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval);
 
     static void
-    integrate(const unsigned int n_components,
-              const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-              VectorizedArrayType *         values_array,
-              VectorizedArrayType *         values_quad,
-              VectorizedArrayType *         gradients_quad,
-              VectorizedArrayType *         hessians_quad,
-              VectorizedArrayType *         scratch_data,
-              const bool                    integrate_values,
-              const bool                    integrate_gradients,
-              const bool                    integrate_hessians,
-              const unsigned int            face_no,
-              const unsigned int            subface_index,
-              const unsigned int            face_orientation,
-              const Table<2, unsigned int> &orientation_map);
+    integrate(
+      const unsigned int                     n_components,
+      const EvaluationFlags::EvaluationFlags integration_flag,
+      VectorizedArrayType *                  values_dofs,
+      FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval);
 
-    static bool
+    static void
     gather_evaluate(
       const unsigned int                          n_components,
-      const std::size_t                           n_face_orientations,
+      const EvaluationFlags::EvaluationFlags      evaluation_flag,
       const Number *                              src_ptr,
       const std::vector<ArrayView<const Number>> *sm_ptr,
-      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-      const MatrixFreeFunctions::DoFInfo &                       dof_info,
-      VectorizedArrayType *                                      values_quad,
-      VectorizedArrayType *                                      gradients_quad,
-      VectorizedArrayType *                                      hessians_quad,
-      VectorizedArrayType *                                      scratch_data,
-      const bool         evaluate_values,
-      const bool         evaluate_gradients,
-      const bool         evaluate_hessians,
-      const unsigned int active_fe_index,
-      const unsigned int first_selected_component,
-      const std::array<unsigned int, VectorizedArrayType::size()> cells,
-      const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
-      const unsigned int                                          subface_index,
-      const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
-      const std::array<unsigned int, VectorizedArrayType::size()>
-                                    face_orientations,
-      const Table<2, unsigned int> &orientation_map);
+      FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval);
 
-    static bool
+    static void
     integrate_scatter(
       const unsigned int                          n_components,
-      const std::size_t                           n_face_orientations,
+      const EvaluationFlags::EvaluationFlags      integration_flag,
       Number *                                    dst_ptr,
       const std::vector<ArrayView<const Number>> *sm_ptr,
-      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-      const MatrixFreeFunctions::DoFInfo &                       dof_info,
-      VectorizedArrayType *                                      values_array,
-      VectorizedArrayType *                                      values_quad,
-      VectorizedArrayType *                                      gradients_quad,
-      VectorizedArrayType *                                      hessians_quad,
-      VectorizedArrayType *                                      scratch_data,
-      const bool         integrate_values,
-      const bool         integrate_gradients,
-      const bool         integrate_hessians,
-      const unsigned int active_fe_index,
-      const unsigned int first_selected_component,
-      const std::array<unsigned int, VectorizedArrayType::size()> cells,
-      const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
-      const unsigned int                                          subface_index,
-      const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
-      const std::array<unsigned int, VectorizedArrayType::size()>
-                                    face_orientations,
-      const Table<2, unsigned int> &orientation_map);
+      FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval);
 
     static bool
     fast_evaluation_supported(const unsigned int given_degree,
index 3f6272b615d21e1f40eb55065c5d4354370c6850..8dabe9598f6090717a7418d6bbcd330a52cef96f 100644 (file)
@@ -44,28 +44,20 @@ namespace internal
   template <int dim, typename Number, typename VectorizedArrayType>
   void
   FEEvaluationFactory<dim, Number, VectorizedArrayType>::evaluate(
-    const unsigned int                                         n_components,
-    const EvaluationFlags::EvaluationFlags                     evaluation_flag,
-    const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
-    VectorizedArrayType *values_dofs_actual,
-    VectorizedArrayType *values_quad,
-    VectorizedArrayType *gradients_quad,
-    VectorizedArrayType *hessians_quad,
-    VectorizedArrayType *scratch_data)
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags evaluation_flag,
+    const VectorizedArrayType *            values_dofs,
+    FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> &eval)
   {
     instantiation_helper_run<
       1,
       FEEvaluationImplEvaluateSelector<dim, VectorizedArrayType>>(
-      shape_info.data[0].fe_degree,
-      shape_info.data[0].n_q_points_1d,
+      eval.get_shape_info().data[0].fe_degree,
+      eval.get_shape_info().data[0].n_q_points_1d,
       n_components,
       evaluation_flag,
-      shape_info,
-      values_dofs_actual,
-      values_quad,
-      gradients_quad,
-      hessians_quad,
-      scratch_data);
+      values_dofs,
+      eval);
   }
 
 
@@ -73,29 +65,21 @@ namespace internal
   template <int dim, typename Number, typename VectorizedArrayType>
   void
   FEEvaluationFactory<dim, Number, VectorizedArrayType>::integrate(
-    const unsigned int                                         n_components,
-    const EvaluationFlags::EvaluationFlags                     integration_flag,
-    const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
-    VectorizedArrayType *values_dofs_actual,
-    VectorizedArrayType *values_quad,
-    VectorizedArrayType *gradients_quad,
-    VectorizedArrayType *hessians_quad,
-    VectorizedArrayType *scratch_data,
-    const bool           sum_into_values_array)
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags integration_flag,
+    VectorizedArrayType *                  values_dofs,
+    FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> &eval,
+    const bool sum_into_values_array)
   {
     instantiation_helper_run<
       1,
       FEEvaluationImplIntegrateSelector<dim, VectorizedArrayType>>(
-      shape_info.data[0].fe_degree,
-      shape_info.data[0].n_q_points_1d,
+      eval.get_shape_info().data[0].fe_degree,
+      eval.get_shape_info().data[0].n_q_points_1d,
       n_components,
       integration_flag,
-      shape_info,
-      values_dofs_actual,
-      values_quad,
-      gradients_quad,
-      hessians_quad,
-      scratch_data,
+      values_dofs,
+      eval,
       sum_into_values_array);
   }
 
@@ -116,40 +100,20 @@ namespace internal
   template <int dim, typename Number, typename VectorizedArrayType>
   void
   FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::evaluate(
-    const unsigned int                                         n_components,
-    const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-    const VectorizedArrayType *                                values_array,
-    VectorizedArrayType *                                      values_quad,
-    VectorizedArrayType *                                      gradients_quad,
-    VectorizedArrayType *                                      hessians_quad,
-    VectorizedArrayType *                                      scratch_data,
-    const bool                                                 evaluate_values,
-    const bool                    evaluate_gradients,
-    const bool                    evaluate_hessians,
-    const unsigned int            face_no,
-    const unsigned int            subface_index,
-    const unsigned int            face_orientation,
-    const Table<2, unsigned int> &orientation_map)
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags evaluation_flag,
+    const VectorizedArrayType *            values_dofs,
+    FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval)
   {
     instantiation_helper_run<
       1,
       FEFaceEvaluationImplEvaluateSelector<dim, VectorizedArrayType>>(
-      data.data[0].fe_degree,
-      data.data[0].n_q_points_1d,
+      eval.get_shape_info().data[0].fe_degree,
+      eval.get_shape_info().data[0].n_q_points_1d,
       n_components,
-      data,
-      values_array,
-      values_quad,
-      gradients_quad,
-      hessians_quad,
-      scratch_data,
-      evaluate_values,
-      evaluate_gradients,
-      evaluate_hessians,
-      face_no,
-      subface_index,
-      face_orientation,
-      orientation_map);
+      evaluation_flag,
+      values_dofs,
+      eval);
   }
 
 
@@ -157,158 +121,70 @@ namespace internal
   template <int dim, typename Number, typename VectorizedArrayType>
   void
   FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::integrate(
-    const unsigned int                                         n_components,
-    const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-    VectorizedArrayType *                                      values_array,
-    VectorizedArrayType *                                      values_quad,
-    VectorizedArrayType *                                      gradients_quad,
-    VectorizedArrayType *                                      hessians_quad,
-    VectorizedArrayType *                                      scratch_data,
-    const bool                                                 integrate_values,
-    const bool                    integrate_gradients,
-    const bool                    integrate_hessians,
-    const unsigned int            face_no,
-    const unsigned int            subface_index,
-    const unsigned int            face_orientation,
-    const Table<2, unsigned int> &orientation_map)
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags integration_flag,
+    VectorizedArrayType *                  values_dofs,
+    FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval)
   {
     instantiation_helper_run<
       1,
       FEFaceEvaluationImplIntegrateSelector<dim, VectorizedArrayType>>(
-      data.data[0].fe_degree,
-      data.data[0].n_q_points_1d,
+      eval.get_shape_info().data[0].fe_degree,
+      eval.get_shape_info().data[0].n_q_points_1d,
       n_components,
-      data,
-      values_array,
-      values_quad,
-      gradients_quad,
-      hessians_quad,
-      scratch_data,
-      integrate_values,
-      integrate_gradients,
-      integrate_hessians,
-      face_no,
-      subface_index,
-      face_orientation,
-      orientation_map);
+      integration_flag,
+      values_dofs,
+      eval);
   }
 
 
 
   template <int dim, typename Number, typename VectorizedArrayType>
-  bool
+  void
   FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::gather_evaluate(
     const unsigned int                          n_components,
-    const std::size_t                           n_face_orientations,
+    const EvaluationFlags::EvaluationFlags      evaluation_flag,
     const Number *                              src_ptr,
     const std::vector<ArrayView<const Number>> *sm_ptr,
-    const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-    const MatrixFreeFunctions::DoFInfo &                       dof_info,
-    VectorizedArrayType *                                      values_quad,
-    VectorizedArrayType *                                      gradients_quad,
-    VectorizedArrayType *                                      hessians_quad,
-    VectorizedArrayType *                                      scratch_data,
-    const bool                                                 evaluate_values,
-    const bool         evaluate_gradients,
-    const bool         evaluate_hessians,
-    const unsigned int active_fe_index,
-    const unsigned int first_selected_component,
-    const std::array<unsigned int, VectorizedArrayType::size()> cells,
-    const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
-    const unsigned int                                          subface_index,
-    const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
-    const std::array<unsigned int, VectorizedArrayType::size()>
-                                  face_orientations,
-    const Table<2, unsigned int> &orientation_map)
+    FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval)
   {
-    return instantiation_helper_run<
+    instantiation_helper_run<
       1,
       FEFaceEvaluationImplGatherEvaluateSelector<dim,
                                                  Number,
                                                  VectorizedArrayType>>(
-      data.data[0].fe_degree,
-      data.data[0].n_q_points_1d,
+      eval.get_shape_info().data[0].fe_degree,
+      eval.get_shape_info().data[0].n_q_points_1d,
       n_components,
-      n_face_orientations,
+      evaluation_flag,
       src_ptr,
       sm_ptr,
-      data,
-      dof_info,
-      values_quad,
-      gradients_quad,
-      hessians_quad,
-      scratch_data,
-      evaluate_values,
-      evaluate_gradients,
-      evaluate_hessians,
-      active_fe_index,
-      first_selected_component,
-      cells,
-      face_nos,
-      subface_index,
-      dof_access_index,
-      face_orientations,
-      orientation_map);
+      eval);
   }
 
 
 
   template <int dim, typename Number, typename VectorizedArrayType>
-  bool
+  void
   FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::integrate_scatter(
     const unsigned int                          n_components,
-    const std::size_t                           n_face_orientations,
+    const EvaluationFlags::EvaluationFlags      integration_flag,
     Number *                                    dst_ptr,
     const std::vector<ArrayView<const Number>> *sm_ptr,
-    const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
-    const MatrixFreeFunctions::DoFInfo &                       dof_info,
-    VectorizedArrayType *                                      values_array,
-    VectorizedArrayType *                                      values_quad,
-    VectorizedArrayType *                                      gradients_quad,
-    VectorizedArrayType *                                      hessians_quad,
-    VectorizedArrayType *                                      scratch_data,
-    const bool                                                 integrate_values,
-    const bool         integrate_gradients,
-    const bool         integrate_hessians,
-    const unsigned int active_fe_index,
-    const unsigned int first_selected_component,
-    const std::array<unsigned int, VectorizedArrayType::size()> cells,
-    const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
-    const unsigned int                                          subface_index,
-    const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
-    const std::array<unsigned int, VectorizedArrayType::size()>
-                                  face_orientations,
-    const Table<2, unsigned int> &orientation_map)
+    FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &eval)
   {
-    return instantiation_helper_run<
+    instantiation_helper_run<
       1,
       FEFaceEvaluationImplIntegrateScatterSelector<dim,
                                                    Number,
                                                    VectorizedArrayType>>(
-      data.data[0].fe_degree,
-      data.data[0].n_q_points_1d,
+      eval.get_shape_info().data[0].fe_degree,
+      eval.get_shape_info().data[0].n_q_points_1d,
       n_components,
-      n_face_orientations,
+      integration_flag,
       dst_ptr,
       sm_ptr,
-      data,
-      dof_info,
-      values_array,
-      values_quad,
-      gradients_quad,
-      hessians_quad,
-      scratch_data,
-      integrate_values,
-      integrate_gradients,
-      integrate_hessians,
-      active_fe_index,
-      first_selected_component,
-      cells,
-      face_nos,
-      subface_index,
-      dof_access_index,
-      face_orientations,
-      orientation_map);
+      eval);
   }
 
 
index b171e4344b2f79fa3a8595e4cd9b4eee4632f5ae..81f6aaa3d4ee13591792b7162a2523f485611cf4 100644 (file)
@@ -2309,7 +2309,6 @@ public:
   void
   integrate(const EvaluationFlags::EvaluationFlags integration_flag);
 
-
   /**
    * @deprecated Please use the integrate() function with the EvaluationFlags argument.
    */
@@ -2329,7 +2328,8 @@ public:
    */
   void
   integrate(const EvaluationFlags::EvaluationFlags integration_flag,
-            VectorizedArrayType *                  values_array);
+            VectorizedArrayType *                  values_array,
+            const bool                             sum_into_values = false);
 
   /**
    * @deprecated Please use the integrate() function with the EvaluationFlags argument.
@@ -2858,6 +2858,8 @@ namespace internal
   inline std::tuple<
     const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *,
     const internal::MatrixFreeFunctions::DoFInfo *,
+    const internal::MatrixFreeFunctions::
+      MappingInfoStorage<dim - is_face, dim, Number, VectorizedArrayType> *,
     unsigned int,
     unsigned int>
   get_shape_info_and_indices(
@@ -2895,6 +2897,7 @@ namespace internal
         active_fe_index,
         active_quad_index),
       &dof_info,
+      &mapping_storage,
       active_fe_index,
       active_quad_index);
   }
@@ -2932,10 +2935,6 @@ inline FEEvaluationBase<dim,
                                                     n_q_points,
                                                     active_fe_index,
                                                     active_quad_index),
-      internal::MatrixFreeFunctions::
-        MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
-          data_in.get_mapping_info(),
-          quad_no),
       quad_no,
       is_interior_face,
       face_type)
@@ -4992,8 +4991,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             for (unsigned int e = 0; e < dim; ++e)
               hessian_out[comp][d][d] +=
                 jac_grad[d][e] *
-                this->gradients_from_hessians_quad[(comp * dim + e) * nqp +
-                                                   q_point];
+                this->gradients_quad[(comp * dim + e) * nqp + q_point];
 
           // add off-diagonal part of J' * grad(u)
           for (unsigned int d = 0, count = dim; d < dim; ++d)
@@ -5001,8 +4999,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
               for (unsigned int f = 0; f < dim; ++f)
                 hessian_out[comp][d][e] +=
                   jac_grad[count][f] *
-                  this->gradients_from_hessians_quad[(comp * dim + f) * nqp +
-                                                     q_point];
+                  this->gradients_quad[(comp * dim + f) * nqp + q_point];
 
           // take symmetric part
           for (unsigned int d = 0; d < dim; ++d)
@@ -5099,8 +5096,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             for (unsigned int e = 0; e < dim; ++e)
               hessian_out[comp][d] +=
                 jac_grad[d][e] *
-                this->gradients_from_hessians_quad[(comp * dim + e) * nqp +
-                                                   q_point];
+                this->gradients_quad[(comp * dim + e) * nqp + q_point];
         }
     }
   return hessian_out;
@@ -7301,28 +7297,13 @@ FEEvaluation<dim,
 
   if (fe_degree > -1)
     SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
-      evaluate(n_components,
-               evaluation_flag_actual,
-               *this->data,
-               const_cast<VectorizedArrayType *>(values_array),
-               this->values_quad,
-               this->gradients_quad,
-               this->hessians_quad,
-               this->scratch_data.begin());
+      evaluate(n_components, evaluation_flag_actual, values_array, *this);
   else
     internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::evaluate(
       n_components,
       evaluation_flag_actual,
-      *this->data,
       const_cast<VectorizedArrayType *>(values_array),
-      this->values_quad,
-      this->gradients_quad,
-      this->hessians_quad,
-      this->scratch_data.begin());
-
-  if (hessians_on_general_cells)
-    for (unsigned int i = 0; i < n_components * dim * n_q_points; ++i)
-      this->gradients_from_hessians_quad[i] = this->gradients_quad[i];
+      *this);
 
 #  ifdef DEBUG
   if (evaluation_flag_actual & EvaluationFlags::values)
@@ -7452,12 +7433,10 @@ namespace internal
   /**
    * Implementation for vectors that have the begin() methods.
    */
-  template <int dim,
-            int fe_degree,
-            int n_q_points_1d,
-            typename Number,
+  template <typename Number,
             typename VectorizedArrayType,
             typename VectorType,
+            typename T,
             typename std::enable_if<
               internal::has_begin<VectorType>::value &&
                 std::is_same<decltype(std::declval<VectorType>().begin()),
@@ -7465,24 +7444,18 @@ namespace internal
               VectorType>::type * = nullptr>
   bool
   try_integrate_scatter_inplace(
+    T                                             phi,
     VectorType &                                  destination,
     const unsigned int                            cell,
-    const unsigned int                            n_components,
     const unsigned int                            active_fe_index,
     const unsigned int                            first_selected_component,
     const internal::MatrixFreeFunctions::DoFInfo *dof_info,
-    VectorizedArrayType *                         values_quad,
-    VectorizedArrayType *                         gradients_quad,
-    VectorizedArrayType *                         hessians_quad,
-    VectorizedArrayType *                         scratch_data,
-    const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *data,
-    const EvaluationFlags::EvaluationFlags integration_flag)
+    const EvaluationFlags::EvaluationFlags        evaluation_flag)
   {
     // If the index storage is interleaved and contiguous and the vector storage
     // has the correct alignment, we can directly pass the pointer into the
-    // vector to the integrate() call, without writing temporary results into a
-    // separate data field that will later be added into the vector. This saves
-    // some operations.
+    // vector to the evaluate() call, without reading the vector entries into a
+    // separate data field. This saves some operations.
     if (std::is_same<typename VectorType::value_type, Number>::value &&
         dof_info->index_storage_variants
             [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] ==
@@ -7505,28 +7478,8 @@ namespace internal
             dof_info->component_dof_indices_offset[active_fe_index]
                                                   [first_selected_component] *
               VectorizedArrayType::size());
-        if (fe_degree > -1)
-          SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
-            integrate(n_components,
-                      integration_flag,
-                      *data,
-                      vec_values,
-                      values_quad,
-                      gradients_quad,
-                      hessians_quad,
-                      scratch_data,
-                      true);
-        else
-          FEEvaluationFactory<dim, Number, VectorizedArrayType>::integrate(
-            n_components,
-            integration_flag,
-            *data,
-            vec_values,
-            values_quad,
-            gradients_quad,
-            hessians_quad,
-            scratch_data,
-            true);
+
+        phi->integrate(evaluation_flag, vec_values, true);
 
         return true;
       }
@@ -7537,31 +7490,23 @@ namespace internal
   /**
    * Implementation for all other vectors like block vectors.
    */
-  template <int dim,
-            int fe_degree,
-            int n_q_points_1d,
-            typename Number,
+  template <typename Number,
             typename VectorizedArrayType,
             typename VectorType,
+            typename T,
             typename std::enable_if<
               !internal::has_begin<VectorType>::value ||
                 !std::is_same<decltype(std::declval<VectorType>().begin()),
                               Number *>::value,
               VectorType>::type * = nullptr>
   bool
-  try_integrate_scatter_inplace(
-    VectorType &,
-    const unsigned int,
-    const unsigned int,
-    const unsigned int,
-    const unsigned int,
-    const internal::MatrixFreeFunctions::DoFInfo *,
-    const VectorizedArrayType *,
-    const VectorizedArrayType *,
-    const VectorizedArrayType *,
-    const VectorizedArrayType *,
-    const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *,
-    const EvaluationFlags::EvaluationFlags)
+  try_integrate_scatter_inplace(T,
+                                VectorType &,
+                                const unsigned int,
+                                const unsigned int,
+                                const unsigned int,
+                                const internal::MatrixFreeFunctions::DoFInfo *,
+                                const EvaluationFlags::EvaluationFlags)
   {
     return false;
   }
@@ -7689,7 +7634,8 @@ FEEvaluation<dim,
              Number,
              VectorizedArrayType>::
   integrate(const EvaluationFlags::EvaluationFlags integration_flag,
-            VectorizedArrayType *                  values_array)
+            VectorizedArrayType *                  values_array,
+            const bool                             sum_into_values_array)
 {
 #  ifdef DEBUG
   if (integration_flag & EvaluationFlags::values)
@@ -7734,24 +7680,16 @@ FEEvaluation<dim,
     SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
       integrate(n_components,
                 integration_flag_actual,
-                *this->data,
                 values_array,
-                this->values_quad,
-                this->gradients_quad,
-                this->hessians_quad,
-                this->scratch_data.begin(),
-                false);
+                *this,
+                sum_into_values_array);
   else
     internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::integrate(
       n_components,
       integration_flag_actual,
-      *this->data,
       values_array,
-      this->values_quad,
-      this->gradients_quad,
-      this->hessians_quad,
-      this->scratch_data.begin(),
-      false);
+      *this,
+      sum_into_values_array);
 
 #  ifdef DEBUG
   this->dof_values_initialized = true;
@@ -7805,22 +7743,13 @@ FEEvaluation<dim,
   integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
                     VectorType &                           destination)
 {
-  if (internal::try_integrate_scatter_inplace<dim,
-                                              fe_degree,
-                                              n_q_points_1d,
-                                              Number,
-                                              VectorizedArrayType>(
+  if (internal::try_integrate_scatter_inplace<Number, VectorizedArrayType>(
+        this,
         destination,
         this->cell,
-        n_components,
         this->active_fe_index,
         this->first_selected_component,
         this->dof_info,
-        this->values_quad,
-        this->gradients_quad,
-        this->hessians_quad,
-        this->scratch_data.begin(),
-        this->data,
         integration_flag) == false)
     {
       integrate(integration_flag, this->begin_dof_values());
@@ -7930,9 +7859,7 @@ FEFaceEvaluation<dim,
       internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
       internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior;
   Assert(this->mapping_data != nullptr, ExcNotInitialized());
-  const internal::MatrixFreeFunctions::FaceToCellTopology<
-    VectorizedArrayType::size()> &faces =
-    this->matrix_info->get_face_info(face_index);
+
   if (face_index >=
         this->matrix_info->get_task_info().face_partition_data.back() &&
       face_index <
@@ -7943,21 +7870,7 @@ FEFaceEvaluation<dim,
              "boundary faces use FEFaceEvaluation with the parameter "
              "is_interior_face set to true. "));
 
-  this->face_no =
-    (this->is_interior_face ? faces.interior_face_no : faces.exterior_face_no);
-  this->subface_index = this->is_interior_face == true ?
-                          GeometryInfo<dim>::max_children_per_cell :
-                          faces.subface_index;
-
-  // First check if interior or exterior cell has non-standard orientation
-  // (i.e. the third bit is one or not). Then set zero if this cell has
-  // standard-orientation else copy the first three bits
-  // (which is equivalent to modulo 8). See also the documentation of
-  // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation.
-  this->face_orientation =
-    (this->is_interior_face == (faces.face_orientation >= 8)) ?
-      (faces.face_orientation % 8) :
-      0;
+  this->reinit_face(this->matrix_info->get_face_info(face_index));
 
   this->cell_type = this->matrix_info->get_mapping_info().face_type[face_index];
   const unsigned int offsets =
@@ -8142,16 +8055,11 @@ FEFaceEvaluation<dim,
   evaluate(const VectorizedArrayType *            values_array,
            const EvaluationFlags::EvaluationFlags evaluation_flag)
 {
-  Assert(
-    (evaluation_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
-                         EvaluationFlags::hessians)) == 0,
-    ExcMessage(
-      "Only EvaluationFlags::values, EvaluationFlags::gradients, and EvaluationFlags::hessians are supported."));
-
-  if (!(evaluation_flag & EvaluationFlags::values) &&
-      !(evaluation_flag & EvaluationFlags::gradients) &&
-      !(evaluation_flag & EvaluationFlags::hessians))
-    return;
+  Assert((evaluation_flag &
+          ~(EvaluationFlags::values | EvaluationFlags::gradients |
+            EvaluationFlags::hessians)) == 0,
+         ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
+                    "and EvaluationFlags::hessians are supported."));
 
   const bool hessians_on_general_cells =
     evaluation_flag & EvaluationFlags::hessians &&
@@ -8160,86 +8068,15 @@ FEFaceEvaluation<dim,
   if (hessians_on_general_cells)
     evaluation_flag_actual |= EvaluationFlags::gradients;
 
-  if (this->dof_access_index ==
-        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-      this->is_interior_face == false)
-    {
-      const auto face_nos          = this->compute_face_no_data();
-      const auto face_orientations = this->compute_face_orientations();
-
-#  ifdef DEBUG
-      // currently on structured meshes are supported -> face numbers and
-      // orientations have to be the same for all filled lanes
-      for (unsigned int v = 1; v < VectorizedArrayType::size(); ++v)
-        {
-          if (face_nos[v] != numbers::invalid_unsigned_int)
-            AssertDimension(face_nos[0], face_nos[v]);
-          if (face_orientations[v] != numbers::invalid_unsigned_int)
-            AssertDimension(face_orientations[0], face_orientations[v]);
-        }
-#  endif
-
-      internal::FEFaceEvaluationImplEvaluateSelector<dim, VectorizedArrayType>::
-        template run<fe_degree, n_q_points_1d>(
-          n_components,
-          *this->data,
-          values_array,
-          this->begin_values(),
-          this->begin_gradients(),
-          this->begin_hessians(),
-          this->scratch_data.begin(),
-          evaluation_flag_actual & EvaluationFlags::values,
-          evaluation_flag_actual & EvaluationFlags::gradients,
-          evaluation_flag_actual & EvaluationFlags::hessians,
-          face_nos[0],
-          this->subface_index,
-          face_orientations[0],
-          this->descriptor->face_orientations);
-    }
+  if (fe_degree > -1)
+    internal::FEFaceEvaluationImplEvaluateSelector<dim, VectorizedArrayType>::
+      template run<fe_degree, n_q_points_1d>(n_components,
+                                             evaluation_flag_actual,
+                                             values_array,
+                                             *this);
   else
-    {
-      if (fe_degree > -1)
-        internal::FEFaceEvaluationImplEvaluateSelector<dim,
-                                                       VectorizedArrayType>::
-          template run<fe_degree, n_q_points_1d>(
-            n_components,
-            *this->data,
-            values_array,
-            this->begin_values(),
-            this->begin_gradients(),
-            this->begin_hessians(),
-            this->scratch_data.begin(),
-            evaluation_flag_actual & EvaluationFlags::values,
-            evaluation_flag_actual & EvaluationFlags::gradients,
-            evaluation_flag_actual & EvaluationFlags::hessians,
-            this->face_no,
-            this->subface_index,
-            this->face_orientation,
-            this->descriptor->face_orientations);
-      else
-        internal::FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::
-          evaluate(n_components,
-                   *this->data,
-                   values_array,
-                   this->begin_values(),
-                   this->begin_gradients(),
-                   this->begin_hessians(),
-                   this->scratch_data.begin(),
-                   evaluation_flag_actual & EvaluationFlags::values,
-                   evaluation_flag_actual & EvaluationFlags::gradients,
-                   evaluation_flag_actual & EvaluationFlags::hessians,
-                   this->face_no,
-                   this->subface_index,
-                   this->face_orientation,
-                   this->descriptor->face_orientations);
-    }
-
-  if (hessians_on_general_cells)
-    {
-      unsigned int size = n_components * dim * n_q_points;
-      for (unsigned int i = 0; i < size; ++i)
-        this->gradients_from_hessians_quad[i] = this->gradients_quad[i];
-    }
+    internal::FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::
+      evaluate(n_components, evaluation_flag_actual, values_array, *this);
 
 #  ifdef DEBUG
   if (evaluation_flag_actual & EvaluationFlags::values)
@@ -8341,26 +8178,21 @@ FEFaceEvaluation<dim,
                  n_components,
                  Number,
                  VectorizedArrayType>::
-  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag,
+  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
             VectorizedArrayType *                  values_array)
 {
-  Assert(
-    (evaluation_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
-                         EvaluationFlags::hessians)) == 0,
-    ExcMessage(
-      "Only EvaluationFlags::values, EvaluationFlags::gradients, and EvaluationFlags::hessians are supported."));
-
-  if (!(evaluation_flag & EvaluationFlags::values) &&
-      !(evaluation_flag & EvaluationFlags::gradients) &&
-      !(evaluation_flag & EvaluationFlags::hessians))
-    return;
+  Assert((integration_flag &
+          ~(EvaluationFlags::values | EvaluationFlags::gradients |
+            EvaluationFlags::hessians)) == 0,
+         ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
+                    "and EvaluationFlags::hessians are supported."));
 
-  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
-  if (evaluation_flag & EvaluationFlags::hessians &&
+  EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
+  if (integration_flag & EvaluationFlags::hessians &&
       (this->cell_type > internal::MatrixFreeFunctions::affine))
     {
       unsigned int size = n_components * dim * n_q_points;
-      if (evaluation_flag & EvaluationFlags::gradients)
+      if (integration_flag & EvaluationFlags::gradients)
         {
           for (unsigned int i = 0; i < size; ++i)
             this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
@@ -8369,43 +8201,19 @@ FEFaceEvaluation<dim,
         {
           for (unsigned int i = 0; i < size; ++i)
             this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
-          evaluation_flag_actual |= EvaluationFlags::gradients;
+          integration_flag_actual |= EvaluationFlags::gradients;
         }
     }
 
   if (fe_degree > -1)
     internal::FEFaceEvaluationImplIntegrateSelector<dim, VectorizedArrayType>::
-      template run<fe_degree, n_q_points_1d>(
-        n_components,
-        *this->data,
-        values_array,
-        this->begin_values(),
-        this->begin_gradients(),
-        this->begin_hessians(),
-        this->scratch_data.begin(),
-        evaluation_flag_actual & EvaluationFlags::values,
-        evaluation_flag_actual & EvaluationFlags::gradients,
-        evaluation_flag_actual & EvaluationFlags::hessians,
-        this->face_no,
-        this->subface_index,
-        this->face_orientation,
-        this->descriptor->face_orientations);
+      template run<fe_degree, n_q_points_1d>(n_components,
+                                             integration_flag_actual,
+                                             values_array,
+                                             *this);
   else
     internal::FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::
-      integrate(n_components,
-                *this->data,
-                values_array,
-                this->begin_values(),
-                this->begin_gradients(),
-                this->begin_hessians(),
-                this->scratch_data.begin(),
-                evaluation_flag_actual & EvaluationFlags::values,
-                evaluation_flag_actual & EvaluationFlags::gradients,
-                evaluation_flag_actual & EvaluationFlags::hessians,
-                this->face_no,
-                this->subface_index,
-                this->face_orientation,
-                this->descriptor->face_orientations);
+      integrate(n_components, integration_flag_actual, values_array, *this);
 }
 
 
@@ -8461,119 +8269,44 @@ FEFaceEvaluation<dim,
     ExcMessage(
       "Only EvaluationFlags::values, EvaluationFlags::gradients, and EvaluationFlags::hessians are supported."));
 
-  const auto fu = [&]() {
-    const auto shared_vector_data = internal::get_shared_vector_data(
-      input_vector,
-      this->dof_access_index ==
-        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
-      this->active_fe_index,
-      this->dof_info);
-
-    if (this->dof_access_index ==
-          internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-        this->is_interior_face == false)
-      {
-        const auto cells             = this->get_cell_or_face_ids();
-        const auto face_nos          = this->compute_face_no_data();
-        const auto face_orientations = this->compute_face_orientations();
+  const auto shared_vector_data = internal::get_shared_vector_data(
+    input_vector,
+    this->dof_access_index ==
+      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
+    this->active_fe_index,
+    this->dof_info);
 
-        return internal::FEFaceEvaluationImplGatherEvaluateSelector<
+  if (fe_degree > 0 &&
+      internal::FEFaceEvaluationImplGatherEvaluateSelector<
+        dim,
+        Number,
+        VectorizedArrayType>::supports(evaluation_flag,
+                                       *this->data,
+                                       internal::get_beginning<Number>(
+                                         input_vector),
+                                       this->dof_info->index_storage_variants
+                                         [this->dof_access_index][this->cell]))
+    {
+      if (fe_degree > -1)
+        internal::FEFaceEvaluationImplGatherEvaluateSelector<
           dim,
           Number,
           VectorizedArrayType>::template run<fe_degree,
                                              n_q_points_1d>(
           n_components,
-          VectorizedArrayType::size(),
+          evaluation_flag,
           internal::get_beginning<Number>(input_vector),
           shared_vector_data,
-          *this->data,
-          *this->dof_info,
-          this->begin_values(),
-          this->begin_gradients(),
-          this->begin_hessians(),
-          this->scratch_data.begin(),
-          evaluation_flag & EvaluationFlags::values,
-          evaluation_flag & EvaluationFlags::gradients,
-          evaluation_flag & EvaluationFlags::hessians,
-          this->active_fe_index,
-          this->first_selected_component,
-          cells,
-          face_nos,
-          this->subface_index,
-          this->dof_access_index,
-          face_orientations,
-          this->descriptor->face_orientations);
-      }
-    else
-      {
-        // TODO: this copying should not be necessary once we have introduced
-        // an internal-data structure
-        std::array<unsigned int, VectorizedArrayType::size()> cells_   = {};
-        std::array<unsigned int, VectorizedArrayType::size()> face_no_ = {};
-        std::array<unsigned int, VectorizedArrayType::size()>
-          face_orientation_ = {};
-
-        cells_[0]            = this->cell;
-        face_no_[0]          = this->face_no;
-        face_orientation_[0] = this->face_orientation;
-
-        if (fe_degree > -1)
-          {
-            return internal::FEFaceEvaluationImplGatherEvaluateSelector<
-              dim,
-              Number,
-              VectorizedArrayType>::template run<fe_degree,
-                                                 n_q_points_1d>(
-              n_components,
-              1,
-              internal::get_beginning<Number>(input_vector),
-              shared_vector_data,
-              *this->data,
-              *this->dof_info,
-              this->begin_values(),
-              this->begin_gradients(),
-              this->begin_hessians(),
-              this->scratch_data.begin(),
-              evaluation_flag & EvaluationFlags::values,
-              evaluation_flag & EvaluationFlags::gradients,
-              evaluation_flag & EvaluationFlags::hessians,
-              this->active_fe_index,
-              this->first_selected_component,
-              cells_,
-              face_no_,
-              this->subface_index,
-              this->dof_access_index,
-              face_orientation_,
-              this->descriptor->face_orientations);
-          }
-        else
-          return internal::
-            FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::
-              gather_evaluate(n_components,
-                              1,
-                              internal::get_beginning<Number>(input_vector),
-                              shared_vector_data,
-                              *this->data,
-                              *this->dof_info,
-                              this->begin_values(),
-                              this->begin_gradients(),
-                              this->begin_hessians(),
-                              this->scratch_data.begin(),
-                              evaluation_flag & EvaluationFlags::values,
-                              evaluation_flag & EvaluationFlags::gradients,
-                              evaluation_flag & EvaluationFlags::hessians,
-                              this->active_fe_index,
-                              this->first_selected_component,
-                              cells_,
-                              face_no_,
-                              this->subface_index,
-                              this->dof_access_index,
-                              face_orientation_,
-                              this->descriptor->face_orientations);
-      }
-  };
-
-  if (!fu())
+          *this);
+      else
+        internal::FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::
+          gather_evaluate(n_components,
+                          evaluation_flag,
+                          internal::get_beginning<Number>(input_vector),
+                          shared_vector_data,
+                          *this);
+    }
+  else
     {
       this->read_dof_values(input_vector);
       this->evaluate(evaluation_flag);
@@ -8633,7 +8366,7 @@ FEFaceEvaluation<dim,
                  n_components_,
                  Number,
                  VectorizedArrayType>::
-  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
+  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
                     VectorType &                           destination)
 {
   Assert((this->dof_access_index ==
@@ -8648,81 +8381,40 @@ FEFaceEvaluation<dim,
     this->active_fe_index,
     this->dof_info);
 
-  // TODO: this copying should not be necessary once we have introduced
-  // an internal-data structure
-  std::array<unsigned int, VectorizedArrayType::size()> cells_            = {};
-  std::array<unsigned int, VectorizedArrayType::size()> face_no_          = {};
-  std::array<unsigned int, VectorizedArrayType::size()> face_orientation_ = {};
-
-  cells_[0]            = this->cell;
-  face_no_[0]          = this->face_no;
-  face_orientation_[0] = this->face_orientation;
-
-  if (fe_degree > -1)
+  if (fe_degree > 0 &&
+      internal::FEFaceEvaluationImplGatherEvaluateSelector<
+        dim,
+        Number,
+        VectorizedArrayType>::supports(integration_flag,
+                                       *this->data,
+                                       internal::get_beginning<Number>(
+                                         destination),
+                                       this->dof_info->index_storage_variants
+                                         [this->dof_access_index][this->cell]))
     {
-      if (!internal::FEFaceEvaluationImplIntegrateScatterSelector<
-            dim,
-            Number,
-            VectorizedArrayType>::template run<fe_degree,
-                                               n_q_points_1d>(
-            n_components,
-            1,
-            internal::get_beginning<Number>(destination),
-            shared_vector_data,
-            *this->data,
-            *this->dof_info,
-            this->begin_dof_values(),
-            this->begin_values(),
-            this->begin_gradients(),
-            this->begin_hessians(),
-            this->scratch_data.begin(),
-            evaluation_flag & EvaluationFlags::values,
-            evaluation_flag & EvaluationFlags::gradients,
-            evaluation_flag & EvaluationFlags::hessians,
-            this->active_fe_index,
-            this->first_selected_component,
-            cells_,
-            face_no_,
-            this->subface_index,
-            this->dof_access_index,
-            face_orientation_,
-            this->descriptor->face_orientations))
-        {
-          // if we arrive here, writing into the destination vector did not
-          // succeed because some of the assumptions in integrate_scatter were
-          // not fulfilled (e.g. an element or degree that does not support
-          // direct writing), so we must do it here
-          this->distribute_local_to_global(destination);
-        }
+      if (fe_degree > -1)
+        internal::FEFaceEvaluationImplIntegrateScatterSelector<
+          dim,
+          Number,
+          VectorizedArrayType>::template run<fe_degree,
+                                             n_q_points_1d>(
+          n_components,
+          integration_flag,
+          internal::get_beginning<Number>(destination),
+          shared_vector_data,
+          *this);
+      else
+        internal::FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::
+          integrate_scatter(n_components,
+                            integration_flag,
+                            internal::get_beginning<Number>(destination),
+                            shared_vector_data,
+                            *this);
     }
   else
     {
-      if (!internal::FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::
-            integrate_scatter(n_components,
-                              1,
-                              internal::get_beginning<Number>(destination),
-                              shared_vector_data,
-                              *this->data,
-                              *this->dof_info,
-                              this->begin_dof_values(),
-                              this->begin_values(),
-                              this->begin_gradients(),
-                              this->begin_hessians(),
-                              this->scratch_data.begin(),
-                              evaluation_flag & EvaluationFlags::values,
-                              evaluation_flag & EvaluationFlags::gradients,
-                              evaluation_flag & EvaluationFlags::hessians,
-                              this->active_fe_index,
-                              this->first_selected_component,
-                              cells_,
-                              face_no_,
-                              this->subface_index,
-                              this->dof_access_index,
-                              face_orientation_,
-                              this->descriptor->face_orientations))
-        {
-          this->distribute_local_to_global(destination);
-        }
+      integrate(integration_flag);
+      this->distribute_local_to_global(destination);
     }
 }
 
index ab1d057382322c930eba31cdcdfe2107c62fdf15..bf7457a1e528db44d773cf429ed6f52f66d65e2f 100644 (file)
@@ -125,6 +125,22 @@ class FEEvaluationBaseData
 public:
   static constexpr unsigned int dimension = dim;
 
+  using NumberType = VectorizedArrayType;
+
+  /**
+   * Constructor, taking data stored in MatrixFree. If applied to problems
+   * with more than one quadrature formula selected during construction of
+   * `matrix_free`, `quad_no` allows to select the appropriate formula.
+   */
+  FEEvaluationBaseData(const std::tuple<const ShapeInfoType *,
+                                        const DoFInfo *,
+                                        const MappingInfoStorageType *,
+                                        unsigned int,
+                                        unsigned int> &info,
+                       const unsigned int              quad_no,
+                       const bool                      is_interior_face,
+                       const unsigned int              face_type);
+
   /**
    * Copy constructor.
    */
@@ -136,6 +152,23 @@ public:
   FEEvaluationBaseData &
   operator=(const FEEvaluationBaseData &other);
 
+  /**
+   * Sets the pointers for values, gradients, hessians to the central
+   * scratch_data_array inside the given scratch array, for a given number of
+   * components as provided by one of the derived classes.
+   */
+  void
+  set_data_pointers(AlignedVector<VectorizedArrayType> *scratch_data,
+                    const unsigned int                  n_components);
+
+  /**
+   * Initialize the indices related to the given face description. Used during
+   * the reinit() call of the FEFaceEvaluation class.
+   */
+  void
+  reinit_face(const internal::MatrixFreeFunctions::FaceToCellTopology<
+              VectorizedArrayType::size()> &face);
+
   /**
    * @name 1: Access to geometry data at quadrature points
    */
@@ -371,26 +404,78 @@ public:
   unsigned int
   get_active_quadrature_index() const;
 
-  //@}
+  /**
+   * If is_face is true, this function returns the face number within a cell for
+   * the face this object was initialized to. On cells where the face makes no
+   * sense, an exception is thrown.
+   *
+   * @note This function depends on the internal representation of data, which
+   * is not stable between releases of deal.II, and is hence mostly for
+   * internal use.
+   */
+  unsigned int
+  get_face_no() const;
 
-protected:
   /**
-   * Constructor. Made protected to prevent users from directly using this
-   * class. Takes all data stored in MatrixFree. If applied to problems with
-   * more than one quadrature formula selected during construction of
-   * `matrix_free`, `quad_no` allows to select the appropriate formula.
+   * If is_face is true, this function returns the index of a subface along a
+   * face in case this object was initialized to a face with hanging nodes. On
+   * faces On cells where the face makes no sense, an exception is thrown.
+   *
+   * @note This function depends on the internal representation of data, which
+   * is not stable between releases of deal.II, and is hence mostly for
+   * internal use.
    */
-  FEEvaluationBaseData(const std::tuple<const ShapeInfoType *,
-                                        const DoFInfo *,
-                                        unsigned int,
-                                        unsigned int> &info,
-                       const MappingInfoStorageType &  mapping_data,
-                       const unsigned int              quad_no,
-                       const bool                      is_interior_face,
-                       const unsigned int              face_type);
+  unsigned int
+  get_subface_index() const;
+
+  /**
+   * If is_face is true, this function returns the orientation index within an
+   * array of orientations as returned by get_orientation_map().
+   *
+   * @note This function depends on the internal representation of data, which
+   * is not stable between releases of deal.II, and is hence mostly for
+   * internal use.
+   */
+  unsigned int
+  get_face_orientation() const;
 
   /**
-   * Constructor that comes with reduced functionality and works similar as
+   * If is_face is true, this function returns the map re-ordering from face
+   * quadrature points in standard orientation to a specific orientation.
+   *
+   * @note This function depends on the internal representation of data, which
+   * is not stable between releases of deal.II, and is hence mostly for
+   * internal use.
+   */
+  const Table<2, unsigned int> &
+  get_orientation_map() const;
+
+  /**
+   * Return the current index in the access to compressed indices.
+   *
+   * @note This function depends on the internal representation of data, which
+   * is not stable between releases of deal.II, and is hence mostly for
+   * internal use.
+   */
+  internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex
+  get_dof_access_index() const;
+
+  /**
+   * Return whether this object was set up to work on what is called
+   * "interior" faces in the evaluation context.
+   *
+   * @note This function depends on the internal representation of data, which
+   * is not stable between releases of deal.II, and is hence mostly for
+   * internal use.
+   */
+  bool
+  get_is_interior_face() const;
+
+  //@}
+
+protected:
+  /**
+   * Constructor that comes with reduced functionality and works similarly as
    * FEValues.
    */
   FEEvaluationBaseData(
@@ -398,15 +483,6 @@ protected:
       internal::MatrixFreeFunctions::
         MappingDataOnTheFly<dim, Number, VectorizedArrayType>> &mapping_data);
 
-  /**
-   * Sets the pointers for values, gradients, hessians to the central
-   * scratch_data_array inside the given scratch array, for a given number of
-   * components as provided by one of the derived classes.
-   */
-  void
-  set_data_pointers(AlignedVector<VectorizedArrayType> *scratch_data,
-                    const unsigned int                  n_components);
-
   /**
    * A pointer to the unit cell shape data, i.e., values, gradients and
    * Hessians in 1D at the quadrature points that constitute the tensor
@@ -683,20 +759,20 @@ template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
   FEEvaluationBaseData(const std::tuple<const ShapeInfoType *,
                                         const DoFInfo *,
+                                        const MappingInfoStorageType *,
                                         unsigned int,
                                         unsigned int> &shape_dof_info,
-                       const MappingInfoStorageType &  mapping_data,
                        const unsigned int              quad_no,
                        const bool                      is_interior_face,
                        const unsigned int              face_type)
   : data(std::get<0>(shape_dof_info))
   , dof_info(std::get<1>(shape_dof_info))
-  , mapping_data(&mapping_data)
+  , mapping_data(std::get<2>(shape_dof_info))
   , quad_no(quad_no)
-  , active_fe_index(std::get<2>(shape_dof_info))
-  , active_quad_index(std::get<3>(shape_dof_info))
+  , active_fe_index(std::get<3>(shape_dof_info))
+  , active_quad_index(std::get<4>(shape_dof_info))
   , descriptor(
-      &mapping_data.descriptor
+      &mapping_data->descriptor
          [is_face ?
             (active_quad_index * std::max<unsigned int>(1, dim - 1) +
              (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
@@ -1081,6 +1157,67 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
 
 
 
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_face_no()
+  const
+{
+  return face_no;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_subface_index() const
+{
+  return subface_index;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_face_orientation() const
+{
+  return face_orientation;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline const Table<2, unsigned int> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_orientation_map() const
+{
+  Assert(descriptor != nullptr, ExcNotInitialized());
+  return descriptor->face_orientations;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_dof_access_index() const
+{
+  return dof_access_index;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline bool
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_is_interior_face() const
+{
+  return is_interior_face;
+}
+
+
+
 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline void
 FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
@@ -1122,6 +1259,29 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
 
 
 
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline void
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::reinit_face(
+  const internal::MatrixFreeFunctions::FaceToCellTopology<
+    VectorizedArrayType::size()> &face)
+{
+  face_no = (is_interior_face ? face.interior_face_no : face.exterior_face_no);
+  subface_index = is_interior_face == true ?
+                    GeometryInfo<dim>::max_children_per_cell :
+                    face.subface_index;
+
+  // First check if interior or exterior cell has non-standard orientation
+  // (i.e. the third bit is one or not). Then set zero if this cell has
+  // standard-orientation else copy the first three bits
+  // (which is equivalent to modulo 8). See also the documentation of
+  // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation.
+  face_orientation = (is_interior_face == (face.face_orientation >= 8)) ?
+                       (face.face_orientation % 8) :
+                       0;
+}
+
+
+
 #endif // ifndef DOXYGEN
 
 
index fc29f2a9d4b0ea91869502deda3e3ad96538eb8f..cf562d5740801e0177961c83db51f3c710e69171 100644 (file)
@@ -29,6 +29,7 @@
 #include <deal.II/fe/mapping_q.h>
 
 #include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/fe_evaluation_base_data.h>
 #include <deal.II/matrix_free/mapping_info.h>
 #include <deal.II/matrix_free/mapping_info_storage.templates.h>
 #include <deal.II/matrix_free/util.h>
@@ -1175,13 +1176,13 @@ namespace internal
           shape_info.dofs_per_component_on_cell;
         constexpr unsigned int hess_dim = dim * (dim + 1) / 2;
 
-        AlignedVector<VectorizedDouble> cell_points(dim * n_mapping_points);
-        AlignedVector<VectorizedDouble> cell_quads(dim * n_q_points);
-        AlignedVector<VectorizedDouble> cell_grads(dim * dim * n_q_points);
-        AlignedVector<VectorizedDouble> cell_grad_grads(dim * hess_dim *
-                                                        n_q_points);
-        AlignedVector<VectorizedDouble> scratch_data(
-          dim * (2 * n_q_points + 3 * n_mapping_points));
+        MappingInfoStorage<dim, dim, double, VectorizedDouble> temp_data;
+        temp_data.copy_descriptor(my_data);
+        FEEvaluationBaseData<dim, double, false, VectorizedDouble> eval(
+          {&shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0);
+
+        AlignedVector<VectorizedDouble> evaluation_data;
+        eval.set_data_pointers(&evaluation_data, dim);
 
         for (unsigned int cell = begin_cell; cell < end_cell; ++cell)
           for (unsigned vv = 0; vv < n_lanes; vv += n_lanes_d)
@@ -1195,7 +1196,7 @@ namespace internal
                   vectorized_load_and_transpose(n_mapping_points * dim,
                                                 plain_quadrature_points.data(),
                                                 start_indices,
-                                                cell_points.data());
+                                                eval.begin_dof_values());
 
                   FEEvaluationFactory<dim, double, VectorizedDouble>::evaluate(
                     dim,
@@ -1203,12 +1204,8 @@ namespace internal
                       (update_flags_cells & update_jacobian_grads ?
                          EvaluationFlags::hessians :
                          EvaluationFlags::nothing),
-                    shape_info,
-                    cell_points.data(),
-                    cell_quads.data(),
-                    cell_grads.data(),
-                    cell_grad_grads.data(),
-                    scratch_data.data());
+                    eval.begin_dof_values(),
+                    eval);
                 }
               if (update_flags_cells & update_quadrature_points)
                 {
@@ -1225,9 +1222,10 @@ namespace internal
                   else
                     for (unsigned int d = 0; d < dim; ++d)
                       for (unsigned int q = 0; q < n_q_points; ++q)
-                        store_vectorized_array(cell_quads[q + d * n_q_points],
-                                               vv,
-                                               quadrature_points[q][d]);
+                        store_vectorized_array(
+                          eval.begin_values()[q + d * n_q_points],
+                          vv,
+                          quadrature_points[q][d]);
                 }
 
               const unsigned int n_points =
@@ -1240,7 +1238,9 @@ namespace internal
                     Tensor<2, dim, VectorizedDouble> jac;
                     for (unsigned int d = 0; d < dim; ++d)
                       for (unsigned int e = 0; e < dim; ++e)
-                        jac[d][e] = cell_grads[q + (d * dim + e) * n_q_points];
+                        jac[d][e] =
+                          eval
+                            .begin_gradients()[q + (d * dim + e) * n_q_points];
 
                     // eliminate roundoff errors
                     if (cell_type[cell] == cartesian)
@@ -1288,13 +1288,13 @@ namespace internal
                           {
                             for (unsigned int e = 0; e < dim; ++e)
                               jac_grad[d][e][e] =
-                                cell_grad_grads[q + (d * hess_dim + e) *
-                                                      n_q_points];
+                                eval.begin_hessians()[q + (d * hess_dim + e) *
+                                                            n_q_points];
                             for (unsigned int c = dim, e = 0; e < dim; ++e)
                               for (unsigned int f = e + 1; f < dim; ++f, ++c)
                                 jac_grad[d][e][f] = jac_grad[d][f][e] =
-                                  cell_grad_grads[q + (d * hess_dim + c) *
-                                                        n_q_points];
+                                  eval.begin_hessians()[q + (d * hess_dim + c) *
+                                                              n_q_points];
                             const auto inv_jac_grad =
                               process_jacobian_gradient(inv_jac,
                                                         inv_jac,
@@ -2088,21 +2088,30 @@ namespace internal
           shape_info.dofs_per_component_on_cell;
         constexpr unsigned int hess_dim = dim * (dim + 1) / 2;
 
-        AlignedVector<VectorizedDouble> cell_points(dim * n_mapping_points);
-        AlignedVector<VectorizedDouble> face_quads(dim * n_q_points);
-        AlignedVector<VectorizedDouble> face_grads(dim * dim * n_q_points);
-        AlignedVector<VectorizedDouble> face_grad_grads(dim * hess_dim *
-                                                        n_q_points);
+        MappingInfoStorage<dim - 1, dim, double, VectorizedDouble> temp_data;
+        temp_data.copy_descriptor(my_data);
+        FEEvaluationBaseData<dim, double, true, VectorizedDouble> eval_int(
+          {&shape_info, nullptr, &temp_data, 0, 0}, 0, true, 0);
+        FEEvaluationBaseData<dim, double, true, VectorizedDouble> eval_ext(
+          {&shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0);
 
-        AlignedVector<VectorizedDouble> scratch_data(
-          dim * (2 * n_q_points + 3 * n_mapping_points));
+        AlignedVector<VectorizedDouble> evaluation_data, evaluation_data_ext;
+        eval_int.set_data_pointers(&evaluation_data, dim);
+        eval_ext.set_data_pointers(&evaluation_data_ext, dim);
 
         for (unsigned int face = begin_face; face < end_face; ++face)
           for (unsigned vv = 0; vv < n_lanes; vv += n_lanes_d)
             {
+              internal::MatrixFreeFunctions::FaceToCellTopology<
+                VectorizedDouble::size()>
+                face_double;
+              face_double.interior_face_no = faces[face].interior_face_no;
+              face_double.exterior_face_no = faces[face].exterior_face_no;
+              face_double.face_orientation = faces[face].face_orientation;
+              face_double.subface_index    = faces[face].subface_index;
+
               // load the geometry field for all SIMD lanes
-              unsigned int       start_indices[n_lanes_d];
-              const unsigned int face_no = faces[face].interior_face_no;
+              unsigned int start_indices[n_lanes_d];
               for (unsigned int v = 0; v < n_lanes_d; ++v)
                 if (faces[face].cells_interior[vv + v] !=
                     numbers::invalid_unsigned_int)
@@ -2111,37 +2120,29 @@ namespace internal
                 else
                   start_indices[v] =
                     faces[face].cells_interior[0] * n_mapping_points * dim;
+
+              eval_int.reinit_face(face_double);
               vectorized_load_and_transpose(n_mapping_points * dim,
                                             plain_quadrature_points.data(),
                                             start_indices,
-                                            cell_points.data());
+                                            eval_int.begin_dof_values());
 
               // now let the matrix-free evaluators provide us with the
               // data on faces
               FEFaceEvaluationFactory<dim, double, VectorizedDouble>::evaluate(
                 dim,
-                shape_info,
-                cell_points.data(),
-                face_quads.data(),
-                face_grads.data(),
-                face_grad_grads.data(),
-                scratch_data.data(),
-                true,
-                true,
-                update_flags_faces & update_jacobian_grads,
-                face_no,
-                GeometryInfo<dim>::max_children_per_cell,
-                faces[face].face_orientation > 8 ?
-                  faces[face].face_orientation - 8 :
-                  0,
-                my_data.descriptor[0].face_orientations);
-
+                EvaluationFlags::values | EvaluationFlags::gradients |
+                  (update_flags_faces & update_jacobian_grads ?
+                     EvaluationFlags::hessians :
+                     EvaluationFlags::nothing),
+                eval_int.begin_dof_values(),
+                eval_int);
 
               if (update_flags_faces & update_quadrature_points)
                 for (unsigned int q = 0; q < n_q_points; ++q)
                   for (unsigned int d = 0; d < dim; ++d)
                     store_vectorized_array(
-                      face_quads[d * n_q_points + q],
+                      eval_int.begin_values()[d * n_q_points + q],
                       vv,
                       my_data.quadrature_points
                         [my_data.quadrature_point_offsets[face] + q][d]);
@@ -2155,7 +2156,9 @@ namespace internal
                 [&](unsigned int                     face_no,
                     int                              is_exterior,
                     unsigned int                     q,
-                    Tensor<2, dim, VectorizedDouble> inv_jac) {
+                    Tensor<2, dim, VectorizedDouble> inv_jac,
+                    FEEvaluationBaseData<dim, double, true, VectorizedDouble>
+                      &eval) {
                   Tensor<2, dim, VectorizedDouble> inv_transp_jac_permut;
                   for (unsigned int d = 0; d < dim; ++d)
                     for (unsigned int e = 0; e < dim; ++e)
@@ -2168,7 +2171,8 @@ namespace internal
                   Tensor<2, dim, VectorizedDouble> jacobi;
                   for (unsigned int e = 0; e < dim; ++e)
                     for (unsigned int d = 0; d < dim; ++d)
-                      jacobi[d][e] = face_grads[(d * dim + e) * n_q_points + q];
+                      jacobi[d][e] =
+                        eval.begin_gradients()[(d * dim + e) * n_q_points + q];
                   Tensor<2, dim, VectorizedDouble> inv_transp_jac =
                     transpose(invert(jacobi));
                   Tensor<3, dim, VectorizedDouble> jac_grad;
@@ -2176,12 +2180,13 @@ namespace internal
                     {
                       for (unsigned int e = 0; e < dim; ++e)
                         jac_grad[d][e][e] =
-                          face_grad_grads[q + (d * hess_dim + e) * n_q_points];
+                          eval.begin_gradients()[q + (d * hess_dim + e) *
+                                                       n_q_points];
                       for (unsigned int c = dim, e = 0; e < dim; ++e)
                         for (unsigned int f = e + 1; f < dim; ++f, ++c)
                           jac_grad[d][e][f] = jac_grad[d][f][e] =
-                            face_grad_grads[q +
-                                            (d * hess_dim + c) * n_q_points];
+                            eval.begin_gradients()[q + (d * hess_dim + c) *
+                                                         n_q_points];
                       const auto inv_jac_grad =
                         process_jacobian_gradient(inv_transp_jac_permut,
                                                   inv_transp_jac,
@@ -2203,6 +2208,7 @@ namespace internal
                 face_type[face] <= affine ? 1 : n_q_points;
               for (unsigned int q = 0; q < n_points_compute; ++q)
                 {
+                  const unsigned int face_no = faces[face].interior_face_no;
                   Tensor<2, dim, VectorizedDouble> jac;
                   for (unsigned int e = 0; e < dim; ++e)
                     {
@@ -2210,7 +2216,9 @@ namespace internal
                         ExtractFaceHelper::reorder_face_derivative_indices<dim>(
                           face_no, e);
                       for (unsigned int d = 0; d < dim; ++d)
-                        jac[d][ee] = face_grads[(d * dim + e) * n_q_points + q];
+                        jac[d][ee] =
+                          eval_int
+                            .begin_gradients()[(d * dim + e) * n_q_points + q];
                     }
                   Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac);
                   for (unsigned int e = 0; e < dim; ++e)
@@ -2227,7 +2235,7 @@ namespace internal
 
                   if (update_flags_faces & update_jacobian_grads)
                     {
-                      compute_jacobian_grad(face_no, 0, q, inv_jac);
+                      compute_jacobian_grad(face_no, 0, q, inv_jac, eval_int);
                     }
 
                   std::array<Tensor<1, dim, VectorizedDouble>, dim - 1>
@@ -2285,47 +2293,43 @@ namespace internal
                       start_indices[v] =
                         faces[face].cells_exterior[0] * n_mapping_points * dim;
 
+                  eval_ext.reinit_face(face_double);
                   vectorized_load_and_transpose(n_mapping_points * dim,
                                                 plain_quadrature_points.data(),
                                                 start_indices,
-                                                cell_points.data());
+                                                eval_ext.begin_dof_values());
 
                   FEFaceEvaluationFactory<dim, double, VectorizedDouble>::
                     evaluate(dim,
-                             shape_info,
-                             cell_points.data(),
-                             face_quads.data(),
-                             face_grads.data(),
-                             face_grad_grads.data(),
-                             scratch_data.data(),
-                             false,
-                             true,
-                             update_flags_faces & update_jacobian_grads,
-                             faces[face].exterior_face_no,
-                             faces[face].subface_index,
-                             faces[face].face_orientation < 8 ?
-                               faces[face].face_orientation :
-                               0,
-                             my_data.descriptor[0].face_orientations);
+                             EvaluationFlags::values |
+                               EvaluationFlags::gradients |
+                               (update_flags_faces & update_jacobian_grads ?
+                                  EvaluationFlags::hessians :
+                                  EvaluationFlags::nothing),
+                             eval_ext.begin_dof_values(),
+                             eval_ext);
 
                   for (unsigned int q = 0; q < n_points_compute; ++q)
                     {
+                      const unsigned int face_no = faces[face].exterior_face_no;
                       Tensor<2, dim, VectorizedDouble> jac;
                       for (unsigned int e = 0; e < dim; ++e)
                         {
                           const unsigned int ee =
                             ExtractFaceHelper::reorder_face_derivative_indices<
-                              dim>(faces[face].exterior_face_no, e);
+                              dim>(face_no, e);
                           for (unsigned int d = 0; d < dim; ++d)
                             jac[d][ee] =
-                              face_grads[(d * dim + e) * n_q_points + q];
+                              eval_ext
+                                .begin_gradients()[(d * dim + e) * n_q_points +
+                                                   q];
                         }
                       Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac);
                       for (unsigned int e = 0; e < dim; ++e)
                         {
                           const unsigned int ee =
                             ExtractFaceHelper::reorder_face_derivative_indices<
-                              dim>(faces[face].exterior_face_no, e);
+                              dim>(face_no, e);
                           for (unsigned int d = 0; d < dim; ++d)
                             store_vectorized_array(
                               inv_jac[ee][d],
@@ -2335,10 +2339,8 @@ namespace internal
 
                       if (update_flags_faces & update_jacobian_grads)
                         {
-                          compute_jacobian_grad(faces[face].exterior_face_no,
-                                                1,
-                                                q,
-                                                inv_jac);
+                          compute_jacobian_grad(
+                            face_no, 1, q, inv_jac, eval_ext);
                         }
 
                       my_data.normals_times_jacobians[1][offset + q] =
index d121e96a06c4d68622a2020d11af091378d95875..4457046663049b1b95b739432d63961346e76f05 100644 (file)
@@ -304,6 +304,17 @@ namespace internal
       unsigned int
       quad_index_from_n_q_points(const unsigned int n_q_points) const;
 
+      /**
+       * Get the descriptor from another instance of QuadratureDescriptor,
+       * clearing all other data fields.
+       */
+      template <typename Number2, typename VectorizedArrayType2>
+      void
+      copy_descriptor(const MappingInfoStorage<structdim,
+                                               spacedim,
+                                               Number2,
+                                               VectorizedArrayType2> &other);
+
       /**
        * Prints a detailed summary of memory consumption in the different
        * structures of this class to the given output stream.
index 69509d10243d2d50d183cf030dc11fb89fbf75bd..aff20a206e964d4d1675cca697ac8bf92cf3b032 100644 (file)
@@ -178,6 +178,47 @@ namespace internal
 
 
 
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    template <typename Number2, typename VectorizedArrayType2>
+    void
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      copy_descriptor(const MappingInfoStorage<structdim,
+                                               spacedim,
+                                               Number2,
+                                               VectorizedArrayType2> &other)
+    {
+      clear_data_fields();
+      descriptor.clear();
+      descriptor.resize(other.descriptor.size());
+      for (unsigned int i = 0; i < descriptor.size(); ++i)
+        {
+          descriptor[i].n_q_points    = other.descriptor[i].n_q_points;
+          descriptor[i].quadrature_1d = other.descriptor[i].quadrature_1d;
+          descriptor[i].quadrature    = other.descriptor[i].quadrature;
+          for (unsigned int d = 0; d < structdim; ++d)
+            {
+              descriptor[i].tensor_quadrature_weights[d].resize(
+                other.descriptor[i].tensor_quadrature_weights[d].size());
+              std::copy(
+                other.descriptor[i].tensor_quadrature_weights[d].begin(),
+                other.descriptor[i].tensor_quadrature_weights[d].end(),
+                descriptor[i].tensor_quadrature_weights[d].begin());
+            }
+          descriptor[i].quadrature_weights.resize(
+            other.descriptor[i].quadrature_weights.size());
+          std::copy(other.descriptor[i].quadrature_weights.begin(),
+                    other.descriptor[i].quadrature_weights.end(),
+                    descriptor[i].quadrature_weights.begin());
+          descriptor[i].face_orientations =
+            other.descriptor[i].face_orientations;
+        }
+    }
+
+
+
     template <int structdim,
               int spacedim,
               typename Number,

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.